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SIMULATION  OF  INFRARED  SPECTRA 

INTRODUCTION 

The  Naval  Research  Laboratory  supports  a  number  of  programs  concerning  the 
infrared  (IR)  spectral  characteristics  of  the  materials  produced  by  the  combustion  of  a 
variety  of  types  of  fuels.  These  programs  include  measurements  of  the  IR  absorption  and 
emission  spectra  of  the  combustion  plumes  of  various  materials  under  differing  conditions. 
In  an  effort  to  more  fully  understand  the  combustion  processes,  a  spectral  simulation 
effort  has  been  initiated.  The  purpose  of  this  effort  is  to  develop  the  capability  of 
simulating  IR  spectra  for  comparison  with  experimental  results.  The  molecules  of  interest 
in  the  simulation  effort  are  atmospheric  constituents  and  their  isotopic  variants. 

Programs  have  been  developed  to  simulate  the  absorption  and  emission  spectra  of  these 
materials  at  different  temperatures  and  under  varying  conditions.  This  report  discusses 
the  structure  of  the  programs  and  the  input  data  which  are  required.  Results  of 
calculations  for  a  variety  of  materials  and  conditions  are  also  presented. 

PROGRAM  STRUCTURE 

Programs  have  been  developed  to  calculate  the  spectra  of  isotopic  variants  of  water 
(H,0,  HDO,  and  D,0),  CO,,  and  HC£  and  DC L  All  programs  are  written  in  FORTRAN 
and  run  on  the  IBM  AT  computer.  In  the  sections  below,  the  theoretical  underpinning  of 
each  set  of  programs  is  discussed.  A  general  discussion  of  the  program  structure  follows 
and  then  the  details  of  program  implementation  are  outlined.  Copies  of  the  program 
codes  are  included  in  the  Appendix.  Also  in  the  Appendix  are  sample  parameter  files 
showing  the  structure  and  formatting  of  the  input  files.  For  both  water  and  CO,,  actual 
input  and  output  files  are  included.  The  discussion  for  each  molecule  assumes  that  an 
absorption  spectrum  is  to  be  calculated.  The  changes  in  the  program  necessary  to 

calculate  an  emission  spectrum  are  detailed. 

Calculations  For  Isotonic  H.,Q 

The  most  difficult  calculation  is  that  for  water  which  is  an  asymmetric  rotor,  king. 
Hainer  and  Cross  originally  determined  expressions  for  the  energv  levels  of  an 


asymmetric  rotor  [I],  Each  energy  level  is  characterized  I"  J.  the  total  angular 
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momentum,  and  by  two  pseudo-quantum  numbers.  K.:  and  Kr  The  pseudo-quantum 
numbers  are  used  to  distinguish  between  the  2J+1  levels  that  have  the  same  value  of  J. 

A  level  is  labelled  by  the  values  of  K  to  which  it  corresponds  in  the  prolate  symmetric 
rotor  (K.:)  and  in  the  oblate  symmetric  rotor  (Kj).  These  two  values  are  combined  to 
define  an  asymmetry  parameter  r 
r  =  K-.j  -  K.j 

The  parameter  r  assumes  the  values  -J,  -J+l,  ....  0,  ....  J-l,  J.  The  energy  levels  are 
ranked  in  order  by  the  value  of  r  with  the  lowest  energy  level  for  a  given  J  being  that 
with  t  *  -J  and  the  highest  that  with  r  =  J.  In  addition,  r  serves  as  an  index  to  the 
parity  of  the  levels.  For  molecules  belonging  to  the  C2v  point  group,  those  levels  with 
even  |r|  are  symmetric  with  respect  to  interchange  of  identical  nuclei  and  those  with 
odd  |  r  |  are  antisymmetric  with  respect  to  interchange  of  identical  nuclei. 

Cross,  Hainer,  and  King  also  derived  a  means  of  calculating  line  strengths  for  the 
lines  in  an  asymmetric  rotor  spectrum  [2],  The  line  strength  calculation  requires  the 
evaluation  of  the  matrix  elements  of  the  direction  cosines  which  are  the  elements  of  an 
orthogonal  transformation  from  molecule-fixed  to  space-fixed  coordinates.  Application  of 
their  method  allows  calculation  of  the  line  strength  of  a  particular  transition  as  a 
function  of  the  asymmetry  of  the  rotor. 

The  starting  point  for  the  programs  described  herein  was  a  program  coded  by  Y. 
Endo  and  kindly  provided  by  J.E.  Butler  (Code  6170.  Chemistry  Division.  Naval  Research 
Laboratory,  Washington,  D.C.).  The  program  was  designed  for  analysis  of  microwave 
spectra  of  asymmetric  rotors.  The  calculation  of  the  transition  frequencies  was  not 
altered,  but  the  intensity  calculation  required  considerable  modification  (e  g.,  in  the 
original  program  the  temperature  could  not  be  varied).  The  dependence  of  the  line 
intensity  S  for  an  absorption  at  the  frequency  u  is  given  by  the  following  expression 
(note  that,  since  only  relative  intensities  from  the  calculation  are  used,  proportionalities 


may  be  used  instead  of  equalities): 

Sa(i  )[l-exp(-c,y  T0)]  g"F  QR(T0)exp(-c2ER"  T0)£2  S°<T0)  i 

v0  *  band  center 
c2  =  second  radiation  constant 
=.  1.438786  cm  K. 

Ta  =  reference  temperature  =  296K. 
g"  =  statistical  weight  of  lower  state 
F  =  vibration  -  rotation  interaction  constant 
Qr(T0)  =  rotational  partition  function  at  T0 

M 

Er=  rotational  energy  of  lower  state 

i2  »  direction  cosines  matrix  element  connecting  lower  and  upper  states 
S°(T0)  =  band  strength  at  T0 
The  band  strength  depends  on  the  following  factors: 

S°  a  V°  lR|  exp  (-c2Gv"/T0)  C 

v  Qv(T0) 

R  =  matrix  element  of  the  rotationless  electric  dipole  moment 
Gv''  =  vibrational  energy  of  lower  state 
QV(T0)  ■  vibrational  partition  function  at  TQ 
Some  comment  is  necessary  on  several  of  the  terms  included  above.  The  term  in  square 
brackets  in  equation  (1)  is  the  induced  emission  term.  It  is  negligibly  small  for  the  case 
in  which  c2i/  >  >  T  but  is  included  for  completeness.  The  statistical  weight  or  nuclear 
statistics  factor  g"  depends  on  the  symmetry  of  the  lower  state  and  will  be  discussed 
further  below.  The  vibration-rotation  interaction  factor  will  be  assumed  to  be  equal  to 
one  in  most  cases.  The  band  strength  for  the  vibrational  transition  is  available  from 
tabulated  values. 

The  temperature  dependence  of  the  line  strength  arises  through  the  rotational  and 
vibrational  Boltzmann  factors  and  the  rotational  and  vibrational  partition  functions,  as 
well  as  a  small  effect  from  the  induced  emission  term.  Fscludinc  the  induced  emission 
term.  McClatchey  [?|  gives  the  temperature  lependence  as 


c«E"  <T-T0) 

i nr 


(?) 


S(T)  =. 


S(T0)  Qv(T0)Qr(Tj) 

exp 

QvlT)  Qr(T) 


where  E"  is  the  lower  state  energy,  including  both  vibrational  and  rotational  energy 
This  dependence  differs  by  a  factor  of  T0,  T  from  that  given  in  some  sources  ([4],  for 

example!  because  the  units  to  be  used  for  S°  are  cm/ molecule  rather  than  cm'2  atm  _1. 

v 

The  temperature  dependence  for  the  bands  in  the  water  molecule  can  be  written  in 
somewhat  simpler  terms  for  two  reasons:  (1)  the  partition  functions  at  T0  are  constants 
and  can  be  omitted  since  only  relative  strengths  are  being  calculated;  and  (2)  the 
vibrational  Boltzmann  factor  is  equal  to  one  because  all  bands  considered  have  the 
vibrational  ground  state  as  the  lower  state.  Note  that  the  temperature  dependence  of 
the  induced  emission  term  can  be  accounted  for  by  multiplying  by  the  ratio  of  that  term 
at  temperature  T  to  that  at  temperature  T0,  cancelling  out  the  T0  term.  The  resulting 
expression  for  the  line  intensity  at  a  temperature  T  is 

Set  ~  [l-exp(-c,t//T)]  - ^ - exp  (-c,ER"/T)  S°(TJ  ‘4| 

The  program  which  calculates  spectra  for  isotopic  water  reads  input  data  from  a 
parameter  file.  This  file  identifies  the  molecule  for  which  the  spectrum  is  to  be 
calculated  and  supplies  the  information  necessary  for  the  simulation;  for  a  complete 
description  of  the  parameter  file,  see  below.  The  program  calculates  the  transition 
frequencies  and  the  line  intensities.  All  lines  with  frequencies  within  the  specified  range 
and  intensities  greater  than  the  chosen  minimum  are  stored.  After  all  lines  have  been 
calculated,  the  stored  lines  are  sorted  in  order  of  increasing  frequency.  The  lines  are 
then  grouped  into  bins  of  equal  width  with  the  summed  intensity  of  all  lines  within  the 
bin  at  the  center  of  the  bin.  A  Gaussian  lineshape  is  applied  to  the  intensity  to  adjust 
the  resolution  for  comparison  with  experimental  spectra.  An  output  file  is  then  written 
for  plotting  the  spectrum. 

The  parameter  file  contains  the  necessars  inputs  for  the  calculation.  The 
parameters  are  discussed  here  in  the  order  in  which  they  appear  in  the  parameter  file 
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See  pages  A2  and  A3  for  sample  parameter  and  input  files.  ISOTOPE  identifies  the 
particular  isotope  under  consideration:  161  for  HoO.  262  for  D.,0.  and  162  for  HDO. 


> 

i 

i 


The  upper  and  lower  state  rotational  constants  and  centrifugal  distortion  constants  are 
read  into  the  array  PAR  (I).  NBINS  is  the  number  of  bins  of  equal  width  into  which  the 
frequency  range  is  divided.  WIDTH  is  the  halfwidth  of  the  Gaussian  lineshnpe  which  is 
applied.  The  variable  IITY  codes  the  band  type  of  the  vibrational  transition  -  100  for 
type  A  transitions,  010  for  type  B  transitions,  and  001  for  type  C  transitions.  The 
minimum  and  maximum  values  of  J  are  JMIN  and  JMAX.  The  maximum  change  allowed  in 
r  is  KTAU.  The  allowed  frequency  range  is  determined  by  the  minimum  and  maximum 
frequency  values  FMIN  and  FMA.X.  The  band  center  of  the  vibrational  transition  is 
FNUZ.  STMIN  is  the  minimum  strength  requirement  for  transitions  to  be  stored.  The 
temperature  for  the  simulation  is  TT  and  the  band  strength  for  the  transition  is  BS.  The 
fundamental  vibrational  frequencies  of  the  molecule  are  FNU1.  FNU2.  and  FNU3.  After 
the  data  have  all  been  read  in,  they  are  printed  to  the  screen  so  that  it  may  be  verified 
that  the  values  are  correct. 


The  program  calculates  frequencies  and  intensities  for  each  branch  in  turn  with  the 
R  branch  first,  followed  by  the  P  branch  and  then  the  Q  branch.  Once  calculated,  the 
frequency  is  tested  to  determine  whether  it  falls  between  the  limits  of  the  frequency 
range.  If  so,  a  level  name  is  determined  and  the  strength  is  calculated.  The  strength 
calculation  is  divided  into  a  number  of  steps.  The  direction  cosines  matrix  element  is 
calculated  and  multiplied  by  the  rotational  Boltzmann  factor  divided  b\  the  rotational 
partition  function.  The  energy  used  in  the  Boltzmann  factor  is  the  initial  state 
rotational  energy  ELI.  The  result  of  this  calculation  is  multiplied  by  the  nuclear 
statistics  factor  to  give  the  parameter  STRENGTH.  The  value  of  the  statistical  weight 
factor  is  determined  by  symmetry  of  the  levels  under  interchange  of  identical  nuclei  and 
is  indexed  by  the  oddness  or  evenness  of  jr|.  Table  I  gives  the  statistical  vveieht  factor 
for  the  isotopes  of  water  As  STRENGTH  for  each  level  is  determined,  the  values  are 
summed  to  give  TOTIN'T.  To  compare  intensities  crrectlv  for  different  isotopes,  the 


intensities  must  be  normalized  by  TOTINT.  This  is  accomplished  by  dividing  the  strength 
of  each  line  by  TOTINT  in  the  SORT  subroutine.  The  calculation  of  the  line  intensity  is 
completed  by  multiplying  STRENGTH  by  the  frequency,  the  induced  emission  term  and  the 
band  strength  then  dividing  by  the  frequency  of  the  band  center  and  the  vibrational 
partition  function.  Note  that  dividing  by  the  band  center  is  required  because  the  band 
center  is  included  in  the  tabulated  S°(T0)  value. 

After  the  line  intensity  is  calculated,  that  value  is  compared  to  STM  IN;  if  it 
exceeds  the  minimum  strength  criterion,  the  frequency,  strength,  and  level  name  are 
stored.  This  process  is  continued  until  all  lines  in  all  three  branches  have  been 
calculated.  The  line  counter  NLINE  is  incremented  each  time  a  line  is  stored.  After  all 
lines  have  been  calculated,  the  lines  are  arranged  in  order  of  increasing  frequency  by  the 
SORT  subroutine.  The  SORT  subroutine  also  performs  the  normalization  by  TOTINT 
described  above.  Because  this  normalization  occurs  after  the  comparison  of  the 
calculated  line  intensity  to  STMIN,  the  minimum  strength  value  in  the  lines  listed  in  the 
program  output  may  not  correspond  directly  to  STMIN.  The  PRINT  subroutine  then 
prints  out  information  for  each  stored  line:  the  band  type;  the  branch;  the  final  and 
initial  values  of  J,  K.1?  and  Kt;  the  transition  frequency;  the  calculated  strength;  and  the 
relative  strength  calculated  by  dividing  each  line  strength  by  the  maximum  line  strength. 


Next  the  subroutine  BINSORT  sorts  ail  lines  into  bins  of  equal  width,  that  width  being 
(FMAX-FMIN),/NBINS.  The  SHAPE  subroutine  applies  a  Gaussian  lineshape  of  halfwidth 
WIDTH  to  the  summed  intensity  at  the  center  of  each  bin.  The  total  intensity  of  each 


bin  is  calculated  as  SMOOTH(I)  after  application  of  the  lineshape.  The  SHAPE  subroutine 
also  writes  an  output  file  containing  FMIN.  FVIA.X.  NBIN'S.  and  SMOOTHU).  1=1  to 
NBINS. 

The  program  described  above  must  be  modified  to  calculate  an  emission  spectrum 
rather  than  an  absorption  spectrum.  There  is  of  course  no  change  in  the  transition 
frequencies.  The  population  of  the  emitting  species  is  that  of  the  upper  state  rather 
than  the  lower  state.  The  energy  used  in  the  rotational  Boltzmann  factor  is  now  the 


iVr. 


final  state  rotational  energy  ELF  rather  than  ELI  (t..e  terminology  of  the  absorption 
process  is  retained;  "final  state"  actually  refers  to  the  upper  state).  Since  the  molecule 
must  also  be  in  the  excited  vibrational  state  to  emit,  a  vibrational  Boltzmann  factor 
exp  t-c,*FNUZ  TT")  must  be  used.  The  induced  emission  term  is  not  used  for  the 
emission  calculation.  The  index  for  the  statistical  weight  factor  is  applied  to  the  upper 
state  rather  than  the  lower  rate.  In  emission,  a  factor  of  i/4  rather  than  u  is  used  in 
the  intensity  calculation  [5].  With  the  exception  of  these  modifications,  the  program  for 
emission  is  the  same  as  that  for  absorption. 

Calculations  For  CO: 

The  calculation  for  a  linear  molecule  is  considerably  simpler  than  that  for  an 
asymmetric  rotor  since  expressions  for  the  energy  and  the  intensity  can  be  expressed  in 
closed  form.  The  expression  for  the  energy  [6]  is 

E(J)=  BJ(J+1)  -  DJ2(J+1)2  +  HJ3(J+n3  i5 

Values  of  the  constants  are  given  in  the  literature.  For  E  states,  the  calculation  is 
straightforward.  In  calculating  energies  for  IT  states,  the  effect  'of  2-type  doubling  must 
be  considered  [6],  Separate  constants  are  tabulated  for  e  and  f  sublevels  of  TI  states  so 
both  manifolds  may  be  calculated. 

In  analyzing  a  particular  transition,  attention  must  be  paid  both  to  which  levels  are 

populated  and  to  the  selection  rules  governing  the  transition.  E+-E+  transitions  are 

u  S 

straightforward  because  the  2-type  doubling  is  absent;  only  P  and  R  branches  appear. 

For  IIU-E*  transitions,  a  Q-branch  appears  in  addition  to  the  P  and  R  branches.  For 
molecules  such  as  C02  with  the  spins  of  the  nuclei  equal  to  zero,  only  the  symmetric 
levels  with  respect  to  interchange  of  equivalent  nuclei  are  occupied.  The  selection  rules 
require  the  combination  of  symmetric  with  symmetric  levels.  The  occupied  levels  in  the 
E  state  are  those  with  even  J  and  the  absorptions  terminate  in  the  e  sublevels  of  the 
upper  state  levels  with  odd  J  for  the  P  and  R  branches.  The  O  branch  transitions 
terminate  in  the  i symmetric)  f  sublevels  of  the  upper  state  levels  with  even  .1.  l  u  a 
n  -n,,  transition.  P.  Q.  and  R  branches  are  allowed.  The  ••  m metric  levels  f-r  i  11  t.ue 
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are  e  sublevels  with  odd  J.  f  sublevels  with  even  J.  Thus  two  sets  ot'  P  and  R  branches 


arise  -  an  e-e  set  with  odd  J  in  the  lower  state  and  an  f-f  set  with  even  J  in  the  lower 
state.  Since  the  Q  branch  transitions  must  connect  two  symmetric  levels,  the  e  <->  f 
selection  rule  is  followed.  The  Q  transitions  with  odd  J  originate  in  e  sublevels  and 
terminate  in  f  sublevels.  The  Q  transitions  with  even  J  originate  in  f  sublevels  and 
terminate  in  e  sublevels. 

The  expression  for  the  intensity  calculation  for  CO.,  bands  is  given  explicitly  bv 


Rothman  [7], 

S  =  S°  Sj£  exp  (-Co  Er"T  0  [l-expl-c,^  T 5 )]  — E - 

0  “  '  Or(T0) 


where 


S°  = 
1/ 


8-3  g“Ia 
3hc  QV(T0) 


;R  - 


exp  i-C2Gv"  T0) 


(  i 


All  terms  are  as  del  .ned  previously;  Sj^  is  the  Honl-London  factor  and  is  the  isotopic 
abundance.  We  will  consider  only  the  isotope  containing  :-C  and  ;<30:  the  constants  will 
be  dropped  since  only  relative  values  are  of  importance.  The  temperature  dependence  is 
introduced  in  the  same  manner  described  above  for  H.,0.  One  complication  in  the  CO., 
case  is  the  fact  that  the  lower  state  is  not  the  ground  v  ibrational  state  in  all  cases. 

That  means  that  the  vibrational  Boltzmann  factor  which  is  included  in  the  band  strength 
may  not  be  equal  to  one.  In  changing  temperatures,  it  is  necessarv  to  divide  by 
exp(-c,Gv"  T0)  and  to  multiply  by  e\p(-c:Gv"  T>  to  calculate  correctly  the  population  in 
the  lower  state  at  a  temperature  T. 

The  initial  section  of  the  program  reads  input  data  from  a  parameter  file.  The 

parameters  are  discussed  here  in  the  order  in  which  they  appear  in  the  parameter  tile 

See  pages  A 40  and  A4I  for  sample  parameter  and  input  files.  The  first  parameter  read 

ID  which  identifies  the  tvpe  of  band  to  be  calculated:  ID=II  for  a  E  -E  band:  10=32 

u  c 

for  a  n,  -E  band;  ID=23  for  a  IT  ~n  band.  If  1 0=  1 1 .  •■nlv  P  and  R  brandies  are 

1  a  -s  '1 

calculated;  if  ID=22  or  ID=23.  P.Q.  and  R  branches  are  calculated  The  next  et  >t 
parameters  contains  the  upper  tnte  constants  Mil  .  1)1  .  and  III  be  n  e.l  ■  i  the 


calculation  ot'  evergy  levels.  For  the  nu-E^  transition,  the  upper  state  energy  lev ei-  are 
those  t'or  the  e  sublevel  t'or  the  P  and  R  branches.  For  the  n?-n.J  transition,  two 
separate  parameter  files  must  be  used  for  P  and  R  branches  -  one  using  the  e  sublevels 
for  odd  J  and  one  using  the  f  sublevels  for  even  J.  The  next  parameters  read  ar^BL  Q. 
DUQ,  and  HUQ,  the  upper  state  levels  used  to  calculate  Q  branch  transitions.  For  the 
nu-Ej  transition,  these  are  for  the  upper  state  f  sublevels.  For  each  of  the  e-e  and  f-f 
pairs  for  the  n  -  ITU  transition,  the  alternate  set  of  sublevels  is  used  for  the  Q  branch 
upper  state  since  the  e<->f  selection  rule  is  followed.  The  final  set  of  parameters  BL. 

DL.  and  HL  is  used  for  the  calculation  of  the  lower  state  energy  levels.  The  next  set  of 
entries  in  the  parameter  file  gives  minimum  and  maximum  values  of  J  iJMIN  and  JMA.Xi 
and  of  the  frequency  (FMIN  and  FMAX).  Note  that  in  the  energy  level  and  transition 
frequency  calculations.  J  is  incremented  by  2.  The  minimum  .alue  of  J  thus  determines 
whether  the  J  values  will  be  odd  or  even.  The  following  three  frequencies  FNU1.  FNL'Z. 
and  FNU3  are  the  fundamental  frequencies  of  the  molecule.  The  next  line  gives  the  band 
center  frequency  (FNUZ)  and  the  vibrational  energies  of  the  lower  (FLOW)  and  upper 
(FHIGHj  states  involved  in  the  transition.  Next  the  band  strength  SV  and  temperature 
TT  are  input.  For  the  n  -  nu  transition,  the  band  strength  is  equally  split  between  the 
e-e  and  f-f  pairs.  The  final  line  in  the  parameter  file  contains  NB1NS  and  WIDTH. 

Following  printing  of  the  input  data  to  the  screen,  three  sets  of  energy  levels  are 
calculated  for  the  indicated  range  of  J:  EUOL  the  upper  state  levels;  ELtl).  the  lower 
state  levels;  and  EUQ(I).  the  upper  state  levels  of  the  other  sublevel  set  for  calculation 
of  the  Q  branch  transitions.  Changes  in  rotational  energies  for  P  and  R  branches 
<EROTP  and  EROTR)  are  calculated;  for  transitions  involving  a  n  state,  the  changes  in 
rotational  energy  for  the  Q  branch  (EROTQ)  are  also  calculated.  The  vibrational 
Boltzmann  factor  and  the  rotational  and  vibrational  partition  functions  are  calculate.!  at 
T,  and  at  T.  For  the  absorption  calculation,  the  vibrational  Boltzmann  factor  uses 
FLOW'  In  a  series  of  loops,  the  frequencies  and  intensities  of  the  individual  liner  in 


there  are  actual  examples  of  input  data  files  HZONL' I  and  COZNL'o  used  for  the  .-.  hand 
of  HoO  and  the  y,  hand  of  CO.,.  The  sample  output  files  generated  h\  running 


ASYMABS  and  C02ABS  respectively  are  also  shown.  These  may  be  useful  as  tests  in 


trying  to  implement  the  programs. 

RESULTS  AND  DISCUSSION 

The  input  parameters  for  the  three  isotopes  of  water  were  collected  from  a  variety 
of  sources,  since  no  one  reference  was  identified  in  which  a  complete  set  for  the  bands 
under  consideration  was  compiled.  The  region  of  the  spectrum  selected  for  consideration 
was  0  -  5000  cm'1.  The  bands  that  make  a  significant  contribution  in  this  region  can  be 
identified  from  the  AFGL  band  strength  listing.  The  most  recent  version  for  the  H20 
bands  was  the  1980  update  [8].  The  four  bands  with  a  band  strength  greater  than  0.5% 
of  the  largest  were  used  in  the  calculation.  These  are  the  three  fundamentals  and  the 
first  overtone  of  the  bending  mode.  The  i/j-k'j  combination  band  has  significant 
intensity  and  is  centered  at  5331  cm'1.  Although  the  lower  frequency  edge  extends 
beyond  5000  cm'1,  this  band  was  not  found  to  make  a  significant  contribution.  The  pure 
rotational  spectrum  was  not  included. 

The  rotational  and  vibrational  parameters  used  for  the  H,0  calculations  are  given  in 
Table  3.  For  the  H,0  calculations,  centrifugal  distortion  constants  are  included  for  all 
bands.  The  corresponding  parameters  for  D;0  are  given  in  Table  4.  but  centrifugal 
distortion  constants  are  not  included.  Note  that  the  band  strengths  are  from  a 
calculation  by  Wilemski  [14]  because  experimental  values  could  not  be  found.  The 
band  is  omitted  because  no  band  strength  could  be  found;  by  analogy  to  HoO  and  HDO.  it 
would  make  a  small  contribution  at  most.  Table  5  contains  the  constants  for  HDO  Note 
that  each  band  has  both  A  and  B  type  character;  separate  band  strengths  are  gi\en  for 
each  contribution. 

The  results  of  calculations  for  absorption  by  H.,0  at  2%K..  I000K.  and  1  "OK.  are 
shown  in  Fig.  1.  The  intensity  is  given  in  arbitrary  units;  the  intensities  of  the 
absorption  spectra  at  different  temperatures  are  correctlv  sealed  relative  to  each  other 
The  width  parameter  here  is  2  cm'1,  corresponding  to  a  resolution  of  4  cm  The 

lowest  frequence  band  at  1594  cm  ;  t-  the  c,  bending  mode  The  very  onall 


contribution  from  2uz  lies  at  3155  cm'1;  it  can  be  seen  in  296K  spectrum.  The  u1  and 
vz  bands  overlie  one  another  in  the  3500-4000  cm'1  region.  As  the  temperature 
increases,  the  spectra  broaden  as  higher  J  levels  are  populated  and  absorb.  The  intensity 
•*of  the  absorption  decreases  as  the  temperature  increases  because  of  the  decrease  in 
population  of  the  absorbing  state.  To  determine  the  accuracy  of  the  calculated  spectra, 
the  frequencies  of  lines  in  the  calculated  spectra  were  compared  to  tabulated  values  for 
each  of  the  bands  studied.  The  calculated  values  agreed  with  literature  values  to  within 
1  cm'1  in  most  cases  with  some  deviations  of  3-4  cm"1  observed.  The  results  were  not 
as  good  for  individual  line  strengths.  Variations  of  50%  or  more  were  noted  in 
comparing  calculated  and  literature  values.  The  purpose  of  this  study,  however,  was  to 
determine  the  region  in  which  emission  occurs  rather  than  the  line-bv-line  strength  of 
particular  portions  of  the  spectrum.  The  accuracy  achieved  is  certainly  adequate  for  that 
purpose. 

Emission  spectra  at  1730K  are  shown  for  all  three  isotopes  in  Fig.  2.  The 
substitution  of  one  D  for  H  in  water  shifts  from  3657  cm'1  to  2724  cm'1,  giving 
three  distinct  bands.  In  going  to  D20,  once  again  u1  and  are  in  the  same  region  but 
all  bands  are  shifted  to  lower  frequency.  Again  all  spectra  are  scaled  correctly  relative 
to  ea  h  other  and  correspond  to  a  resolution  of  4  cm'1.  It  is  clear  from  Fig.  2  that 
changing  from  H,0  to  D,0  shifts  the  most  intense  absorption  band  into  a  significantly 
different  region  of  the  spectrum.  Fig.  3  shows  a  composite  emission  spectrum  at  1730K 
for  an  equimolar  mixture  of  H,0  and  D,0. 

Table  6  gives  the  rotational  and  vibrational  constants  for  C02.  The  bands  included 
are  those  five  for  which  S^>I0-4  cm/ molecule.  Three  of  these  are  S-S  transitions  with 
one  fl-E  transition  and  one  17-17  transition.  The  lower  state  for  the  n -  IT  transition  is  not 
the  ground  vibrational  state  but  instead  is  the  01  101  state  (for  band  notation,  see  [3]). 
Absorption  spectra  calculated  at  296K.  I000K.  and  l~30K  are  shown  in  Fig.  4.  correcth 
scaled  relative  to  one  another.  The  most  prominent  hand  at  2340  cm'1  is  the  i>3 
fundamental.  The  vz  fundamental  appears  at  607  cm'1  and  small  contributions  are  made 
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by  combination  bands  in  the  3600-3700  cm'1  region.  As  the  temperature  is  increased,  the 
absorption  decreases  because  ot'  the  smaller  number  of  molecules  in  the  lower  vibrational 
state.  The  relative  importance  of  the  2300  cm'1  peak  increases  relative  to  the  667  cm'1 
peak  because  of  the  contribution  made  to  the  higher  frequency  peak  by  the  n-n 
transition  originating  in  the  01 101  state.  This  excited  lower  state  has  a  higher 
population  at  the  higher  temperature  and  the  absorption  due  to  this  transition  therefore 
is  increased.  Figure  5  shows  emission  due  to  CO,  at  1000K.  and  1730K.  with  correct 
relative  scaling.  The  2300  cm'1  feature  is  again  the  predominant  one  with  its  importance 
relative  to  the  667  cm'1  band  increasing  with  increasing  temperature.  The  contribution 
due  to  the  3600-3700  cm'1  peaks  is  small  even  at  1730K.  Figure  6  shows  on  an  expanded 
scale  the  main  CO,  band  with  emission  calculated  at  1730K.  and  absorption  calculated  at 
296K..  The  maximum  value  of  each  spectrum  is  individually  scaled  to  a  value  of  10.  If  a 
material  producing  hot  CO,  were  burned  and  the  spectrum  observed  through  an 
atmospheric  path,  the  cold  CO,  in  the  atmosphere  would  absorb  the  radiation  making  up 
the  central  portion  of  the  emission  spectrum.  The  observed  spectrum  would  show  only 
the  wings  of  the  hot  C02  spectrum,  as  shown  in  the  lower  panel  of  Fig.  6. 

Table  7  gives  the  rotational  and  vibrational  constants  used  to  calculate  the  H Ci 
and  DC£  spectra.  The  calculation  also  used  a  vibration-rotation  interaction  constant 
expression  [22]: 

F=  1  -  2.5599xl0'2  m  +  3.203xl0'4  m2  (8) 

Emission  spectra  for  HCf  and  D Ci  at  1730K  are  shown  in  Fig.  7.  Again  it  is  clear  that 
deuteration  of  HC£  shifts  the  emission  into  a  different  region  of  the  spectrum. 

CONCLUSIONS 

Spectral  simulation  programs  have  been  developed  which  correctly  calculate  the 
spectra  of  species  produced  in  the  combustion  of  a  variety  of  fuels.  The  programs 
calculate  both  positions  and  intensities  in  the  v  ibrational- rotational  spectra  of  isotopic 
torms  ot  H,0.  CO,,  and  HCf.  The  calculated  values  of  the  frequencies  agree  well  with 


tabulated  values  while  the  agreement  ot'  intensities  with  literature  values  is  not  as  good. 

As  expected,  the  regions  of  emission  and  absorption  shift  markedly  with  isotopic 
substitution.  The  relative  importance  of  particular  transitions  varies  with  temperature. 

The  simulated  spectra  can  be  combined  to  mimic  observation  through  an  atmospheric  path 
where  both  emission  and  absorption  due  to  the  same  species  are  important. 
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appendix 


The  Appendix  contains  program  listings,  sample  inputs,  and  sample  outputs  for  the  programs 
described  in  the  preceding  report.  Contents  of  the  Appendix  are  listed  below. 


Listing  Page 

Sample  Parameter  File  for  Water  a2 

Sample  Input  File  H20NU1  A.' 

Output  from  ASYMABS  a -l 

Program  ASYMABS  a o 

Program  ASYMEMS  A  I  S 

Subroutines  ASYMSUB  a  30 

Subroutine  SHAPE  A3~ 

Subroutine  BINSORT  AVI 

Sample  Parameter  File  for  C02  A40 

Sample  Input  File  C02NU3  *  \4I 

Output  from  C02ABS  \ 42 

Program  C02ABS  \44 

Program  C02EMS  A4X 


A  5  2 


Program  HCLEMS 


SAMPLE  PARAMETER  FILE  FOR  WATER 


ISOTOPE 

PAR(l)  PAR ( 2 )  PAR ( 3 ) 

PAR (4 )  PAR (5)  PAR (6) 

PAR(7 )  PAR ( 8 )  PAR ( 9 ) 

PAR( 16)  PAR (17)  PAR ( 18 ) 

PAR( 19)  PAR (20)  PAR (21) 

PAR (22)  PAR (23)  PAR(24) 

NBINS [ 15 ]  WIDTH [ F10 . 6 ] 

IITY [ 14 ]  JM-N[I4]  JMAX[I4]  KT2U[I4] 

FMIN [ F10 . 4 ]  FMAX[ F10 . 4 ]  FNUZ[F10.4] 

STMIN [E12 . 3 ]  TT [F10 . 4 ]  BS[F10.4] 

FNU1[F10.4]  FNU2 [ F10 . 4 ]  FNU3[F10.4] 


ISOTOPE  »  isotope  code 

PAR ( 1 )  , PAR ( 2 ) , PAR ( 3 )  ■  A,B,C  constants  for  upper  state 
PAR (4)  -  PAR ( 8 )  ■  distortion  constants  for  upper  state 
PAR (16) , PAR (17) , PAR (18)  *  A,B,C  constants  for  lower  state 
PAR (19)  -  PAR (23)  =  distortion  constants  for  lower  state 
NBINS  *  number  of  bins 
WIDTH  =  Gaussian  halfwidth 
IITY  =  band  type 

JMIN,JMAX  =  minimum  and  maximum  values  of  J 
KTAU  =  maximum  change  in  TAU 

FMIN , FMAX  =  minimum  and  maximum  values  of  frequency 
FNUZ  »  band  center 
STMIN  =  minimum  strength 
TT  =  temperature 
BS  =  band  strength 

FNU1 , FNU2 , FNU3  =  fundamental  frequencies 


SAMPLE  INPUT  FILE  H20NU1 


ft}! 
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1 


TO 


i'll 


:.'*M 


'*9 


m 


§ 


(>Vi 


.»;y 

$s 

K 


J 


$ 

4 


1 

$ 


<K 

cs 

iS 

i 


\Oi 


i 

sf 


m 


9 


& 


161 

27.12217 
0.001233 
0.00049987 
27.880678 
0.00124894 
0.00050838 
5000  2.0 

010  0  2 

0000. 

1.50 

3657.053  1 


14 .3047 

-0.0053874 

0.0012405 

14.521689 

-0.0057655 

0.0013007 


9.1045 

0.03023 

0.0 

9.277459 

0.0325199 

0.0 


25  5 

5000. 
296. 


3657.053 

48.62 


1594.778  3755.93 


.  »v.y  yy 


OUTPUT  FROM  PROGRAM  ASYMABS.FOR 


PARAMETER  FILE  *  H20NU1 


OUTPUT  FILE  3  H20ANU1 


ISOTOPE  *  161 


UPPER  STATE  CONSTANTS 

A  *  27.1222 
DU  3  0.00123300 
DSJ  *  0.00049987 


LOWER  STATE  CONSTANTS 

A  *  27.8807 
DLJ  *  0.00124894 
DSJ  *  0.00050838 


8  *  14.3047 
DLJK*  -0.00538740 
OSK  3  0.00124050 


8  3  14.5217 
OLJK3  -0.00576550 
OSK  =  0.00130070 


C  3  9.1045 

DLK  3  0.03023000 


C  3  9.2775 

DLK  3  0.03251990 


SELECTION  RULES  F  T  F 

J  RANGE  0  TO  25 


MAX.  DELTA  TAU 
FREO.  RANGE 
BAND  ORIGIN 
NU  #1 
NU  #2 
NU  #3 

MIN.  STRENGTH 
TEMPERATURE 
NBINS 
WIDTH 


5 

0.0000  TO  5000.0000  CM- 1 
3657.0530  CM-1 
3657.0530  CM-1 
1594.7780  CM-1 
3755.9300  CM-1 
1.500000 
296.00 
5000 

2.000000 


TOTAL  NUMBER  OF  LINES  = 


a 


TOTAL  BAND  INTENSITY  *  7.5561E-01 


B 

p 

3 

0 

3 

4 

1 

4 

.  3566.55133 

2.288E-00 

0.6224 

B 

p 

2 

1 

2 

3 

0 

3 

3598. 12016 

2 . 048E+00 

0.5572 

B 

p 

1 

0 

1 

2 

1 

2 

3600.96185 

2.330E-00 

0.6340 

B 

a 

3 

1 

2 

3 

2 

1 

3615.21502 

2.070E+00 

0.5632 

B 

a 

3 

0 

3 

3 

1 

2 

3618.02353 

2. 131E»00 

0.5797 

B 

0 

1 

0 

1 

1 

1 

0 

3638.08552 

2.820E+00 

0.7672 

B 

a 

1 

1 

0 

1 

0 

1 

3674.65685 

3. 118E+00 

0.8482 

B 

Q 

3 

1 

2 

3 

0 

3 

3690.60424 

2.634E+00 

0.7165 

B 

a 

3 

2 

1 

3 

1 

2 

3691.35115 

2.591E-00 

0.7049 

a 

R 

2 

1 

2 

1 

0 

1 

3711.08686 

3. 149E+00 

0.8566 

B 

R 

3 

0 

3 

2 

1 

2 

3711.89401 

2.816E*00 

0.7662 

B 

R 

4 

1 

4 

3 

0 

3 

3740.84666 

3.676E+00 

1.0000 

B 

R 

2 

2 

1 

1 

1 

0 

3746.18036 

2.904E+00 

0.7900 

a 

R 

5 

0 

5 

4 

1 

4 

3751.52891 

3. 135E+00 

0.8529 

a 

R 

6 

1 

6 

5 

0 

5 

3770.55376 

2.453E+00 

0.6674 

B 

R 

4 

2 

3 

3 

1 

2 

3777.88815 

2.037E+00 

0.5541 

a 

R 

3 

3 

0 

2 

2 

1 

3800.25831 

2.956E+00 

0.8043 

a 

R 

4 

3 

2 

3 

2 

1 

3818.54410 

2.003E+00 

0.5449 

a 

R 

4 

4 

1 

3 

3 

0 

3846.98903 

2.068E+00 

0.5626 

0,  5000 

TOTAL  BAND  INTENSITY (GAUSS  I  AN  SUM)  *  4.9228E+ 


MAXIMUM  INTENSITY  =  1 .3810E+00 


vV» 


nooooo  ooo^o  nnoo  ooonoooo  non 


PROGRAM 


PROGRAM  ASYMABS . FOR 

IMPLICIT  REAL* 8 ( A-H ,  0-Z ) 

ASYMMETRIC  ROTOR  FREQUENCY  ESTIMATION 
FOR  INFRARED  AND  VISIBLE  SPECTROSCOPY 

CODED  BY  Y.E. 

MODIFIED  BY  C .  H.  D.  ,  R.  T.  L.  ,  H .  H.  N . 


DIMENSION  NAME ( 18 ) ,HD(50) ,HF(50) ,HDLO(50) ,HFLO(50) , 

1  EL(50) , ELLO ( 50 ) ,CFAC(50) ,CFACP(50) ,TVF(50) , 

2  TVI (50), LTYPE ( 3 )  , P ( 3  0 ) 

LOGICAL  LTYPE 

COMMON  /CON/  PAR (30) 

COMMON  /FREQ/  FMIN , FMAX 

COMMON/BINS/  BINST(5000) , BINDEX ( 5000 ) , NBINS , WIDTH 
COMMON/ OTPUT/ACCES 1 $ 

COMMON/ INTEN/TOTINT , STRENGTH 
CHARACTER* 2 4  ACCESS? , ACCES1? 

CHARACTER* 8  FILENM? , FILNM1? 

INTEGER* 2  I YEAR, IMONTH, IDAY , IHOUR, IMINUTE 
CALL  GETTIM ( IHOUR , IMINUTE ) 

CALL  GETDAT ( I YEAR , IMONTH , IDAY ) 


ENTER  NAME  OF  PARAMETER  FILE  TO  BE  USED  AS  INPUT 

WRITE  (6,*)  'ENTER  INPUT  FILE  NAME?' 

READ  99 , FILENM? 

FORMAT (A8) 

ENTER  NAME  OF  FILE  TO  BE  USED  FOR  STORAGE  OF  OUTPUT 

WRITE  (6,*)  'ENTER  OUTPUT  FILE  NAME?' 

READ  99 , FILNM1? 

THE  INPUT  PARAMETER  FILE  IS  LOCATED  ON  DRIVE  C:  IN 
SUBDIRECTORY  ASYMRTR 

THE  OUTPUT  FILE  IS  WRITTEN  TO  DRIVE  A: 

ACCESS?  =  ' C : \ASYMRTR\ ' //FILENM?// ' . DAT ' 

ACCES1?  =  'A: ’//FILNMl?// ' . DAT' 

OPEN  ( 2, FILE=ACCESS?,ACCESS=' SEQUENTIAL' ) 

C 

C  ISOTOPE  CODE  IS  READ: 

C  161  =  H20 

C  262  =  D20 

C  162  =  HDO 

C 
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READ (2,*)  ISOTOPE 


C 

C  THE  VARIABLES  ARE  READ  FROM  THE  PARAMETER  FILE: 

C  PAR (1,2,3)  =  A ,  B ,  C  VALUES  FOR  THE  UPPER  STATE 

C  PAR ( 4 , 5 , 6 , 7 , 8 )  =  CENTRIFUGAL  DISTORTION  CONSTANTS 

C  FOR  THE  UPPER  STATE 

C  PAR (10-15)  =  UNUSED  IN  PRESENT  VERSION;  ALL  =  0 

C  PAR (16, 17, 18)  =  A , B , C  VALUES  FOR  THE  LOWER  STATE 

C  PAR(19, 20, 21, 22, 23)  =  CENTRIFUGAL  DISTORTION 

C  CONSTANTS  FOR  THE  LOWER  STATE 

C  PAR (24-30)  =  UNUSED  IN  PRESENT  VERSION;  ALL  =  0 

C 

DO  111  1=1,7, 3 

111  READ (2 , *)  (PAR(J) ,  J=I,I+2) 

67  FORMAT (3F16. 8) 

DO  121  1=10,15 

121  PAR (I) =0 . 

DO  131  1=16,22,3 

131  READ ( 2 , * )  (PAR(J) ,  J=I,I+2) 

DO  141  1=25,30 
141  PAR ( I ) =0 . 

C 

C  READ  NBINS  =  NUMBER  OF  BINS  INTO  WHICH  FREQUENCY  RANGE  IS 

C  DIVIDED  AND  WIDTH  =  HALF-WIDTH  OF  GAUSSIAN  LINE  SHAPE 

C 

READ  (2,84)  NBINS, WIDTH 
84  FORMAT  (I5,F10.6) 

C 

C  READ  IITY  =  BAND  TYPE 

C  100  =  A  TYPE  BAND 

C  010  =  B  TYPE  BAND 

C  001  =  C  TYPE  BAND 

C  JMIN  AND  JMAX  SET  MINIMUM  AND  MAXIMUM  VALUES  OF  J 

C  KTAU  =  MAXIMUM  CHANGE  ALLOWED  IN  TAU 

C 

READ  (2,69)  IITY , JMIN , JMAX , KTAU 
69  FORMAT  (414) 

C 

C  FMIN  AND  FMAX  SET  ALLOWED  LIMITS  FOR  FREQUENCY  RANGE 

C  IN  CM-1.  FNUZ  IS  THE  BAND  CENTER  IN  CM-1 

C 

READ  (2,79)  FMIN , FMAX , FNUZ 
C 

C  STMIN  =  MINIMUM  STRENGTH  FOR  STORAGE  OF  LINE 

C  TT  =  TEMPERATURE 

C  BS  =  BAND  STRENGTH  FOR  TRANSITION  BEING  CALCULATED 

C 

READ (2, 59)  STMIN, TT , BS 
59  FORMAT (E12.3,2F10.4) 

C 

C  FNU1 , FNU2 , FNU3  ARE  THE  FUNDAMENTAL  FREQUENCIES  IN  CM-1 

C 

READ  (2,79)  FNU1 , FNU2 , FNU3 
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79 


FORMAT (3F10. 4) 
CLOSE  (2) 


THE  PARAMETERS  READ  FROM  THE  PARAMETER  FILE  ARE 
PRINTED  OUT. 


6019 


6022 


6023 


6029 

6034 

6036 

6020 


WRITE  (6,6019) 

FORMAT  (10X,' OUTPUT  FROM  PROGRAM  ASYMABS . FOR '/// ) 

WRITE (6, 6022)  IMONTH, IDAY, I YEAR 
FORMAT (IX, 12,  ' / '  , 12 ,  ' /  '  , 14 ) 

WRITE (6, 6023)  IHOUR, IMINUTE 
FORMAT ( IX , I 2 , ' : ’ ,12/) 

WRITE (6, 6034)  FILENM$ 

WRITE (6,6036)  FILNM1$ 

WRITE (6, 6029)  ISOTOPE 
FORMAT ('  ISOTOPE  =',I4//) 

FORMAT ('  PARAMETER  FILE  =  ' ,A8//)  . 

FORMAT ('  OUTPUT  FILE  =  ’,A8//) 

WRITE(6, 6020)  (PAR(I),  1=1,8) 

FORMAT ('  UPPER  STATE  CONSTANTS'// 

L  '  A  = ' , F16 . 4 , 5X , 1 B  = ' , F16 . 4 , 5X , ' C  =',F16.4/ 

2  '  DLJ  = ' , F16 . 8 , 5X, ' DLJK= ' , F16 . 8 , 5X, ' DLK  =',F16.3/ 

3  '  DSJ  = 1 , F16 . 8 , 5X , ' DSK  =',F16.8//) 

WRITE  (6,6025)  (PAR(I),  1=16,23) 
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6025  FORMAT ( '  LOWER  STATE  CONSTANTS'// 

1  '  A  = ' , F16 . 4 , 5X, ' B  =' ,F16.4,5X, 'C  =',F16.4/ 

2  '  DU  =  '  ,  F16 . 8 , 5X  ,  '  DUK=  '  ,  F16 . 8 , 5X  ,  '  DLK  =',F16.3/ 

3  '  DSJ  = ' , F16 . 8 , 5X, ' DSK  =',F16.8//) 

DO  100  1=1,3 
LTYPE ( I )  = . FALSE . 

IF  (IITY.GE. 100)  LTYPE (1)=. TRUE. 

IF  (MOD(IITY, 100) .GE. 10)  LTYPE ( 2 )=. TRUE . 

IF  (MOD ( IITY , 10 )  .GE.  1)  LTYPE ( 3 )=. TRUE . 

WRITE (6, 6030)  LTYPE , JMIN , JMAX 
FORMAT (//'  SELECTION  RULES  ',3L4/ 

1  '  J  RANGE  ' , 14 , '  TO  ' , 14) 

WRITE (6,6040)  KTAU , FMIN , FMAX , FNUZ , FNU1 , FNU2 , FNU3 , 

1  STMIN, TT,NBINS, WIDTH, BS 


6040  FORMAT (// 


MAX.  DELTA  TAU 
FREQ.  RANGE 
BAND  ORIGIN 
NU  41 
NU  42 
NU  4  3 

MIN.  STRENGTH 

TEMPERATURE 

NBINS 

WIDTH 

'  BAND  STRENGTH 


MV 

'  , F10 . 4 ,  '  TO  ’ 
'  , F10 . 4 ,  '  CM-1 
'  , F10 . 4 ,  '  CM-1 
'  , F10 . 4 ,  ’  CM-1 
'  ,  F10 . 4  ,  '  CM-1 
'  ,  E12 . 3/ 

'  ,  F7 . 2/ 

M5/ 

'  ,  F10 . 6/ 

' , F10 .4/// 


TO  ' , F10 . 4 , ’  CM-1 ' 
CM-1 '  / 

CM-1  '/ 

CM-1  ’/ 

CM-1  '/ 


SET  LINE  COUNTER  TO  ZERO 


NLINE  =  0 


C 

C  SET  INITIAL  TOTAL  BAND  INTENSITY  =  0 

C 

TOTINT=0 . 0 
C 

C  DEFINE  FCT  =  FACTOR  FOR  BOLTZMANN  CALCULATIONS 

C 

FCT=-1 . 438786/TT 
C 

C  CALCULATE  VIBRATIONAL  PARTITION  FUNCTION 

C 

VIBPF-1./ ( ( 1 - -DEXP ( FCT*FNU1) ) * ( 1 . -DEXP ( FCT*FNU2 ) ) * 

1  (1. -DEXP ( FCT* FNU3) ) ) 

C 

C  CALCULATE  ROTATIONAL  PARTITION  FUNCTION 

C 

ROTPF= (TT**3/ ( PAR ( 16 ) *PAR(17) *PAR(18) ) ) **. 5 
C 
C 

DO  9000  J=JMIN , JMAX 
C 

C  PRINT  J  ON  SCREEN  TO  TRACK  PROGRESS  OF  CALCULATION 

C 

PRINT  6010, J 
6010  FORMAT  (IX, 15) 

C 

C 

C******  R  BRANCH  TRANSITIONS 
C 

IF  (J.EQ.JMIN)  GO  TO  4000 
FLJIN=1..  DO/DBLE  ( J) 

JLO  =  J-l 
C 

DO  2000  NTYPE=1 , 3 

IF  ( . NOT • LTYPE ( NTYPE ) )  GO  TO  2000 
C 
C 

DO  2100  ISYMF=1 , 4 
C 

ISYMI  =  ISEL ( ISYMF, NTYPE, 1) 

CALL  ECALC ( 0 , ISYMF, J , HD , HF , EL, NHF , KMAXF , MAXT) 

CALL  ECALC ( 1 , ISYMI , JLO , HDLO , HFLO , ELLO , NHI , KMAXI , MAXTLO ) 
IF  (NHF*NHI . EQ . 0)  GO  TO  2100 
C 

C******  CALCULATE  DIRECTION  COSINE 
C 

CALL  CFACTR ( 1 , NTYPE , NHF , NHI , KMAXF , KMAXI , J , CFAC , CFACP) 

C 

C******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT-MAXTLO 

NSMAX  =  (KTAU-NTDEL+100) 74-25 


NSMIN  =  2  5- (KTAU+NTDEL+100) / 4 
C 
C 

DO  2200  ISHIFT=NSMIN , NSMAX 
IF  (ISHIFT. GE. 0)  GO  TO  2180 
N  =  MINO (NHF+ ISHIFT , NHI ) 

IF  (N.LT.l)  GO  TO  2200 
NPLUSF  =  -ISHIFT 
NPLUSI  =  0 
GO  TO  2210 
C 

2180  N  *  MINO (NHF, NHI- ISHIFT) 

IF  (N.LT.l)  GOTO  2200 
NPLUSF  =  0 
NPLUSI  =  ISHIFT 
2210  CONTINUE 
C 

C******  CALCULATE  INDIVIDUAL  TRANSITIONS 
C 

DO  2300  1=1 , N 
INOF  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL (INOF) 

ELI  =  ELLO (INOI) 

FREQ  =  (ELF-ELI) /I . +FNUZ 
C 

C  TEST  WHETHER  LINE  POSITION  IS  WITHIN  SPECIFIED 

C  FREQUENCY  RANGE.  IF  SO,  CONTINUE  WITH  INTENSITY 

C  CALCULATION;  IF  NOT , CALCULATE  NEXT  LINE. 

C 

IF  ( FREQ. GT.FMAX. OR. FREQ. LT.FMIN)  GO  TO  2300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2* (1-INOF) 

KCF  =  J+l-KAF-ISYMF+2  * ( ISYMF/2 ) 

KAI  =  KMAXI+2* ( 1-INOI) 

KCI  =  JLO+l-KAI-ISYMI+2 * (ISYMI/2) 

C 

C ******  CALCULATE  EIGEN  VECTORS 
C 

CALL  EIVEC ( NHF , HD , HF , ELF , TVF ) 

CALL  EIVEC (NHI , HDLO , HFLO , ELI , TVI ) 

C 

C******  CALCULATE  LINE  STRENGTH 
C 

FAC  =  0.0D0 
IPLUS  =  0 

IF  (KMAXF-KMAXI.EQ. 2)  IPLUS  =  1 
NN  =  MINO (NHF-IPLUS, NHI) 

DO  2520  11=1, NN 

2520  FAC  =  FAC+CFAC ( II ) *TVF ( 11+ I PLUS ) *TVI ( II ) 

IF  (NTYPE.EQ.l)  GOTO  2540 


o o o o  ooooo  ooonoo  non  nonnoon  ooon 


FACP  =  0.0D0 

NN  =  MINO (NHF-1 , NHI ) 

IF  (NN.LT.l)  GO  TO  2512 
DO  2530  11=1, NN 

2530  FACP  =  FACP+CFACP ( II ) *TVF ( II+l ) *TVI ( II ) 

2512  IF  (NTYPE.EQ.2)  FACP  =  -FACP 
FAC  =  FAC+FACP 
2540  STRE  =  FAC**2*FLJIN 

MULTIPLY  DIRECTION  COSINES  FACTOR  BY  ROTATIONAL  BOLTZMANN 
FACTOR  OVER  ROTATIONAL  PARTITION  FUNCTION 

STRE= (STRE*DEXP (-ELI/ (0.695*TT) ) ) /ROTPF 

ASSIGN  NUCLEAR  STATISTICS  FACTOR  BASED  ON  IDENTITY  OF 
ISOTOPE  AND  WHETHER  LEVEL  IS  ODD  OR  EVEN 

H20  1:3  FOR  EVEN: ODD 

D20  6:3  FOR  EVEN: ODD 

HDO  1  FOR  ALL  LEVELS 

KK=IABS (KAI-KCI) 

IF ( ISOTOPE. EQ. 161)  THEN 
STNUC= ( 1 . +2 . * (REAL (MOD (KK, 2 ) ) ) ) 

ELSEIF( ISOTOPE. EQ. 2 62)  THEN 
STNUC= ( 6 • -3 . * (REAL (MOD (KK, 2) ) ) ) 

ELSE 

STNUC=1 . 0 
ENDIF 

STIMULATED  EMISSION  CONTRIBUTION  IS  CALCULATED. 

STIMEM= ( 1 . -DEXP (FCT*FREQ) ) 

CALCULATE  STRENGTH  WHICH  INCLUDES  ROTATIONAL  AND  NUCLEAR 
STATISTICS  FACTORS.  THIS  WILL  BE  USED  LATER  TO  NORMALIZE 
FOR  COMPARISON  OF  DIFFERENT  ISOTOPES.  TOTINT  IS  SUM  OF 
STRENGTH  FOR  EACH  LINE. 

STRENGTH=STNUC*STRE 

CALCULATE  FINAL  INTENSITY  BY  MULTIPLYING  BY  THE  INDUCED 
EMISSION  TERM,  THE  FREQUENCY , AND  THE  BAND  STRENGTH  AND 
DIVIDING  BY  VIBRATIONAL  PARTITION  FUNCTION  AND  BAND  CENTER. 

STRE= (STRENGTH*STIMEM*FREQ*BS ) / ( FNUZ*VIBPF) 

******  STORE  THE  CALCULATED  TRANSITION  IF  STRE  EXCEEDS  THE 
MINIMUM  STRENGTH  CRITERION. 

IF  ( STRE . GE . STMIN)  CALL  STORE ( NLINE , J , KAF , KCF , JLO , KAI , 

1  KCI, FREQ, STRE, NTYPE) 

2300  CONTINUE 
2200  CONTINUE 
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2100  CONTINUE 
2000  CONTINUE 
C 
C 

c******  p  BRANCH  TRANSITIONS  -  P  BRANCH  TRANSITION  FREQUENCIES 
C  AND  STRENGTHS  ARE  CALCULATED  BY  THE  SAME  METHOD  USED 

C  ABOVE  FOR  THE  R  BRANCH  TRANSITIONS. 

C 

4000  CONTINUE 

C 

C 

IF  (J.EQ.JMIN)  GO  TO  4200 
FLJIN  =  1 . DO/DBLE ( J) 

JLO  =  J-l 
C 
C 

DO  3000  NTYPE=1 , 3 

IF  ( . NOT. LTYPE (NTYPE) )  GO  TO  3000 
C 

DO  3100  ISYMF=1 , 4 

ISYMI  =  ISEL (IS YMF, NTYPE, 1) 

CALL  ECALC ( 1 , IS YMF , J , HD , HF , EL , NHF , KMAXF , MAXT ) 

CALL  ECALC ( 0 , ISYMI , JLO , HDLO , HFLO , ELLO , NHI , KMAXI , MAXTLO ) 
IF  (NHF*NHI . EQ . 0 )  GO  TO  3100 
C 

C******  CALCULATE  DIRECTION  COSINES 
C 

CALL  CFACTR ( 1 , NTYPE , NHF , NHI, KMAXF , KMAXI , J , CFAC , CFACP ) 

C 

C******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT -MAXTLO 
NSMAX  =  (KTAU-NTDEL+100) /4-25 
NSMIN  =  25— (KTAU+NTDEL+100) / 4 
C 
C 

DO  3200  ISHIFT=NSMIN, NSMAX 
IF  (ISHIFT. GE.O)  GO  TO  3180 
N  =  MINO (NHF+ISHIFT,NHI) 

IF  (N.LT.l)  GOTO  3200 
NPLUSF  »  -ISHIFT 
NPLUSI  a  o 
GOTO  3210 
C 

3180  N  =  MINO (NHF, NHI-ISHIFT) 

IF  (N.LT.l)  GOTO  3200 
NPLUSF  =  o 
NPLUSI  =  ISHIFT 
3210  CONTINUE 
C 

C******  CALCULATE  INDIVIDUAL  TRANSITIONS 
C 

DO  3300  1=1, N 
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INOr  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL(INOF) 

ELI  =  ELLO(INOI) 

FREQ  =  (ELI -ELF) /I . +FNUZ 

IF  (FREQ.GT. FMAX.OR. FREQ. LT. FMIN)  GO  TO  3300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2* ( 1-INOF) 

KCF  =  J+l-KAF-ISYMF+2* ( ISYMF/2 ) 

KAI  =  KMAXI+2* (1-INOI) 

KCI  =  JLO+l-KAI-ISYMI+2* (ISYMI/2) 

C 

C******  CALCULATE  EIGEN  VECTORS 
C 

CALL  EIVEC (NHF , HD, HF , ELF , TVF) 

CALL  EIVEC ( NHI , HDLO , HFLO , ELI , TVI ) 

C 

c******  CALCULATE  LINE  STRENGTH 
C 

FAC  =  0.0D0 
I PLUS  =  0 

IF  (KMAXF-KMAXI . EQ . 2 )  IPLUS  =  I 
NN  =  MINO (NHF-IPLUS , NHI) 

DO  3520  11=1, NN 

3520  FAC  =  FAC+CFAC ( II ) *TVF ( II+IPLUS ) *TVI ( II ) 

IF  (NTYPE.EQ. 1)  GOTO  3540 

FACP  =  0.0D0 

NN  =  MINO ( NHF- 1, NHI) 

IF  (NN. LT. 1)  GO  TO  3512 
DO  3530  11=1, NN 

3530  FACP  =  FACP+CFACP (II) *TVF (II+l) *TVI ( II ) 

3512  IF  (NTYPE.EQ. 2)  FACP  =  -FACP 
FAC  =  FAC+FACP 
3540  STRE  =  FAC**2*FLJIN 

STRE= (STRE*DEXP ( -ELI/ ( 0 . 695*TT) ) ) /ROTPF 
KK=IABS (KAI-KCI) 

IF ( ISOTOPE. EQ. 161)  THEN 
STNUC= ( 1 . +2 . * ( REAL ( MOD ( KK , 2 ) ) ) ) 

ELSEIF ( ISOTOPE. EQ. 262)  THEN 
STNUC= (6.-3.* ( REAL (MOD (KK, 2 ) ) ) ) 

ELSE 

STNUC=1 . 0 
ENDIF 

STIMEM= ( 1 . -DEXP ( FCT*  FREQ ) ) 

STRENGTH=STNUC*STRE 

STRE= (STRENGTH*STIMEM*FREQ*BS ) / (FNUZ*VIBPF) 

C 

c******  STORE  THE  CALCULATED  TRANSITION 
C 

IF  (STRE.GE.STMIN)  CALL  STORE (NLINE , JLO , KAI , KCI , J , KAF , 
1  KCF, FREQ, STRE, NTYPE) 


c 

3300  CONTINUE 
3200  CONTINUE 
3100  CONTINUE 
3000  CONTINUE 
C 
C 

4100  CONTINUE 
C 

c 

c******  q  BRANCH  TRANSITIONS  -  Q  BRANCH  TRANSITION  FREQUENCIES 
C  AND  LINE  STRENGTHS  ARE  CALCULATED  BY  THE  SAME  METHOD 

C  USED  ABOVE  FOR  THE  P  AND  R  BRANCHES. 

C 

IF  (J.EQ.O)  GO  TO  4200 

FLJIN  =  DBLE (J+J+l) /DBLE ( J* (J+l) ) 

C 

C 

DO  7000  NTYPE=1 , 3 

IF  ( . NOT . LTYPE ( NTYPE ) )  GO  TO  7000 
C 
C 

DO  7100  ISYMF=1 , 4 

ISYMI  =  ISEL(ISYMF, NTYPE, 0) 

CALL  ECALC ( 0 , IS YMF , J , HD , HF , EL , NHF , KMAXF , MAXT ) 

CALL  ECALC ( 1 , ISYMI , J , HDLO , HFLO , ELLO , NHI , KMAXI , MAXTLO ) 

IF  (NHF*NHI . EQ . 0 )  GO  TO  7100 
C 

C******CALCULATE  DIRECTION  COSINES 
C 

CALL  CFACTR ( 0 , NTYPE , NHF , NHI , KMAXF , KMAXI , J , CFAC , CFACP) 

C 

C******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT -MAXTLO 
NSMAX  *  (KTAU-NTDEL+100) / 4-25 
NSMIN  *  25- (KTAU+NTDEL+100) /4 
C 
C 

DO  7200  ISHIFT=NSMIN, NSMAX 
C 

IF  (ISHIFT. GE.O)  GO  TO  7320 
N  =  MINO (NHF+ ISHIFT , NHI ) 

IF  (N.LT. 1)  GO  TO  7200 
NPLUSF  =  -ISHIFT 
NPLUSI  =  0 
GO  TO  7330 
C 

7320  N  =  MINO (NHF, NHI-ISHIFT) 

IF  (N.LT. 1)  GO  TO  7200 
NPLUSF  =  0 
NPLUSI  =  ISHIFT 
C 


c******  INDIVIDUAL  TRANSITIONS 
C 

7330  DO  7300  I=1,N 

INOF  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL (INOF) 

ELI  =  ELLO(INOI)  * 

FREQ  =  ( ELF-ELI ) / 1 . +FNUZ 

IF  ( FREQ • GT . FMAX . OR . FREQ . LT . FMIN )  GO  TO  7300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2* ( 1-INOF) 

KCF  =  J+l-KAF-MOD ( IS YMF , 2 ) 

KAI  =  KMAXI+2  * ( 1-INOI ) 

KCI  =  J+l-KAI-MOD ( ISYMI , 2 ) 

C 

C******EIGEN  vectors 

c 

CALL  EIVEC(NHF,HD,HF,ELF,TVF) 

CALL  EIVEC (NHI , HDLO , HFLO , ELI , TVI ) 

c 

c******line  strength 

c 

FAC  =  0.0D0 

NN  =  MINO (NHI, NHF) 

DO  7520  11=1, NN 


7520 


7514 

C 

7510 


7516 

C 

7512 

C 

7540 


FAC  =  FAC+CFAC (II) *TVF ( II ) *TVI (II) 

IF  (NTYPE.EQ. 1)  GO  TO  7540 

FACP  =  0.0D0 

IF  (KMAXF.LT.J)  GO  TO  7510 
NN  =  MINO ( NHF- 1 , NHI ) 

IF  (NN.LT.l)  GO  TO  7512 
DO  7514  11=1, NN 

FACP  =  FACP+CFACP ( II ) *TVF ( II+l) *TVI ( II ) 
GOTO  7512 

NN  =  MINO (NHF, NHI-1) 

IF  (NN.LT.l)  GO  TO  7512 
DO  7516  11=1, NN 

FACP  =  FACP+CFACP ( II ) *TVF( II) *TVI (II+l) 


IF  (NTYPE.EQ. 3) 
FAC  =  FAC+FACP 


FACP  =  -FACP 


STRE  =  FAC**2*FLJIN 

STRE= (STRE+DEXP ( -ELI/ (0 . 695*TT) ) ) /ROTPF 
KK=IABS (KAI-KCI) 

IF (ISOTOPE . EQ .161)  THEN 
STNUC= ( 1 . +2 . * (REAL (MOD ( KK , 2 ) ) ) ) 
ELSEIF(ISOTOPE. EQ. 262)  THEN 
STNUC= ( 6 . -3 . * ( REAL ( MOD ( KK , 2 ) ) ) ) 


ELSE 

STNUC=1. 0 
ENDIF 

STIMEM= ( 1 . - DEXP ( FCT  * FREQ ) ) 

STRENGTH=STNUC  *  STRE 

STRE= ( STRENGTH*  STIMEM*  FREQ  *  BS ) / ( FNUZ  * VI BPF ) 

C 

C******STORE  THE  CALCULATED  TRANSITION 
C 

IF  (STRE.GE.STMIN)  CALL  STORE (NLINE , J , KAF , KCF , J , KAI , 

1  KCI, FREQ, STRE, NTYPE) 

C 

7300  CONTINUE 
7200  CONTINUE 
7100  CONTINUE 
7000  CONTINUE 
C 
C 

4200  CONTINUE 
C 
C 
C 

9000  CONTINUE 

C 

C 

C******  SORT  AND  PRINTOUT  THE  CALCULATED  SPECTRUM 
C 

C  PRINT  OUT  THE  TOTAL  NUMBER  OF  LINES  MEETING  THE  STRENGTH 

C  AND  FREQUENCY  RANGE  CRITERIA. 

C 

WRITE (6, 6013)  NLINE 

6013  FORMAT (10X, 'TOTAL  NUMBER  OF  LINES  «  ',16/) 

C 

C  ARRANGE  THE  STORED  TRANSITIONS  IN  ORDER  OF  INCREASING 

C  FREQUENCY . 

C 

CALL  SORT (NLINE) 

C 

C  PRINT  OUT  THE  TRANSITION  TYPE  AND  QUANTUM  STATE  LABELS, 

C  THE  TRANSITION  FREQUENCY,  THE  LINE  STRENGTH  AND 

C  THE  RELATIVE  STRENGTH.  DISABLE  THIS  SUBROUTINE  CALL  TO 

C  SUPPRESS  PRINTING. 

C 

C  CALL  PRINT (NLINE) 

C 

C  PRINT  OUT  THE  SUM  OF  STRENGTH  WHICH  INCLUDES  ONLY  THE 

C  DIRECTION  COSINE  AND  THE  ROTATIONAL  TERMS. 

C 


WRITE  (6,6017)  TOTINT 

6017  FORMAT (1 OX, 'TOTAL  BAND  INTENSITY  =  ’,1PE10.4/) 

C 

C  SORT  THE  STORED  TRANSITIONS  INTO  FREQUENCY  BINS  OF  WIDTH 

C  EQUAL  TO  THE  FREQUENCY  RANGE  DIVIDED  BY  NBINS . 


•+ 


3: 


o  o  o  o 


c 

CALL  BINSORT (NLINE) 

APPLY  A  GAUSSIAN  LINE  SHAPE  OF  HALFWIDTH  =  "WIDTH"  TO  THE 
SUMMED  INTENSITY  IN  EACH  BIN. 

CALL  SHAPE (NLINE) 

C 

STOP 

END 
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J 


nno  oooooo  oonvo  ooon  nnoooooo  non 


PROGRAM  ASYMEMS . FOR 
IMPLICIT  REAL* 8 ( A-H , 0-Z ) 

ASYMMETRIC  ROTOR  FREQUENCY  ESTIMATION  PROGRAM 
FOR  INFRARED  AND  VISIBLE  SPECTROSCOPY 

CODED  BY  Y.E. 

MODIFIED  BY  C. H . D. , R. T. L. , H. H. N. 


DIMENSION  NAME ( 18 ) ,HD(50) ,HF(50) ,HDLO(50) ,HFLO(50) , 

1  EL(50) , ELLO ( 50 ) ,CFAC(50) ,CFACP(50) ,TVF(50) , 

2  TVI(50) ,LTYPE(3) ,P(30) 

LOGICAL  LTYPE 

COMMON  /CON/  PAR (30) 

COMMON  /FREQ/  FMIN , FMAX 

COMMON/BINS/  BINST(5000) , BINDEX ( 5000 ) , NBINS , WIDTH 
COMMON/OTPUT/ ACCES I $ 

COMMON/ INTEN/TOTINT , STRENGTH 
CHARACTER* 2 4  ACCESS? , ACCES1$ 

CHARACTER* 8  FILENM? , FILNM1? 

INTEGER* 2  I YEAR , IMONTH , I DAY , IHOUR, IMINUTE 
CALL  GETTIM ( IHOUR , IMINUTE ) 

CALL  GETDAT(IYEAR, IMONTH, IDAY) 


ENTER  NAME  OF  PARAMETER  FILE  TO  BE'  USED  AS  INPUT 

WRITE  (6,*)  'ENTER  INPUT  FILE  NAME?' 

READ  99 , FILENM? 

FORMAT (A8) 

ENTER  NAME  OF  FILE  TO  BE  USED  FOR  STORAGE  OF  OUTPUT 

WRITE  (6,*)  'ENTER  OUTPUT  FILE  NAME?' 

READ  99 , FILNM1? 

THE  INPUT  PARAMETER  FILE  IS  LOCATED  ON  DRIVE  C:  IN 
SUBDIRECTORY  ASYMRTR 

THE  OUTPUT  FILE  IS  WRITTEN  TO  DRIVE  A: 

ACCESS?  =  ' C: \ASYMRTR\ '//FILENM?// ' . DAT ' 

ACCES 1?  =  'A: '//FILNM1?//' .DAT' 

OPEN  (2 ,FILE=ACCESS?,ACCESS=' SEQUENTIAL' ) 

ISOTOPE  CODE  IS  READ: 

161  =  H20 
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C  262  =  D20 

C  162  =  HDO 

C 

READ (2,*)  ISOTOPE 
C 

C  THE  VARIABLES  ARE  READ  FROM  THE  PARAMETER  FILE: 

C  PAR (1,2, 3)  =  A, B, C  VALUES  FOR  THE  UPPER  STATE 

C  PAR (4, 5, 6, 7, 8)  =  CENTRIFUGAL  DISTORTION  CONSTANTS 

C  FOR  THE  UPPER  STATE 

C  PAR (10-15)  =  UNUSED  IN  PRESENT  VERSION;  ALL  =  0 

C  PAR (16, 17, 18)  =  A , B , C  VALUES  FOR  THE  LOWER  STATE 

C  PAR( 19, 2 0,2 1,2 2, 23)  =  CENTRIFUGAL  DISTORTION 

C  CONSTANTS  FOR  THE  LOWER  STATE 

C  PAR (24-30)  =  UNUSED  IN  PRESENT  VERSION;  ALL  =  0 

C 

DO  111  1=1, 7, 3 

111  READ ( 2 , * )  (PAR(J) ,  J=I,I+2) 

67  FORMAT ( 3F16 . 8 ) 

DO  121  1=10,15 

121  PAR ( I ) =0 . 

DO  131  1=16,22,3 

131  READ ( 2 , * )  ( PAR ( J ) ,  J=I,I+2) 

DO  141  1=25,30 
141  PAR ( I ) =0 . 

C 

C  READ  NBINS  =  NUMBER  OF  BINS  INTO  WHICH  FREQUENCY  RANGE  IS 

C  DIVIDED  AND  WIDTH  =  HALF -WIDTH  OF  GAUSSIAN  LINE  SHAPE 

C 

READ  (2,84)  NBINS, WIDTH 
84  FORMAT  (I5,F10.6) 

C 

C  READ  IITY  =  BAND  TYPE 

C  100  =  A  TYPE  BAND 

C  010  =  B  TYPE  BAND 

C  001  =  C  TYPE  BAND 

C  JMIN  AND  JMAX  SET  MINIMUM  AND  MAXIMUM  VALUES  OF  J 

C  KTAU  =  MAXIMUM  CHANGE  ALLOWED  IN  TAU 

C 

READ  (2,69)  IITY , JMIN , JMAX , KTAU 
69  FORMAT  (414) 

C 

C  FMIN  AND  FMAX  SET  ALLOWED  LIMITS  FOR  FREQUENCY  RANGE 

C  IN  CM-1.  FNUZ  IS  THE  BAND  CENTER  IN  CM-1 

C 

READ  (2,79)  FMIN , FMAX , FNUZ 

58  FORMAT ( 3F10 . 6) 

C 

C  STMIN  =  MINIMUM  STRENGTH  FOR  STORAGE  OF  LINE 

C  TT  =  TEMPERATURE 

C  BS  =  BAND  STRENGTH  FOR  TRANSITION  BEING  CALCULATED 

C 

READ (2, 59)  STMIN, TT,BS 

59  FORMAT (E12.3,2F10.4) 
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c 

C  FNU1 , FNU2 , FNU3  ARE  THE  FUNDAMENTAL  FREQUENCIES  IN  CM-1 

C 

READ  (2,79)  FNU1 , FNU2 , FNU3 
79  FORMAT (3F10 . 4 ) 

CLOSE  (2) 

C 

C  THE  PARAMETERS  READ  FROM  THE  PARAMETER  FILE  ARE 

C  PRINTED  OUT. 

C 

WRITE  (6,6019) 

6019  FORMAT  (10X, 'OUTPUT  FROM  PROGRAM  ASYMEMS . FOR ’ /// ) 

WRITE (6, 6022)  IMONTH , IDAY , IYEAR 

6022  F0RMAT(1X,I2, '/• ,12, '/' ,14) 

WRITE (6, 6023)  IHOUR, IMINUTE 

6023  FORMAT (IX, 12 ,  ' :  '  , 12/) 

WRITE (6, 6034)  FILENM$ 

WRITE (6,6036)  FILNM1$ 

WRITE (6, 6029)  ISOTOPE 

6029  FORMAT ('  ISOTOPE  =',I4//) 

6034  FORMAT ( '  PARAMETER  FILE  =  • ,A8//) 

6036  FORMAT ( '  OUTPUT  FILE  =  ' ,A8//) 

WRITE(6, 6020)  (PAR(I),  1=1,8) 

6020  FORMAT ( '  UPPER  STATE  CONSTANTS'// 

1  '  A  = 1 , F16 . 4 , 5X , ' B  =' ,F16.4,5X, 'C  =',F16.4/ 

2  '  DLJ  = ' , F16 . 8 , 5X, ' DLJK= 1 , F16 . 8 , 5X, ' DLK  =',F16.3/ 

3  '  DSJ  = ' , F16 . 8 , 5X , ' DSK  =',F16.8//) 

WRITE  (6,6025)  (PAR(I),  1=16,23) 

6025  FORMAT ('  LOWER  STATE  CONSTANTS'// 

1  '  A  =’ ,F16.4,5X, 'B  «',F16.4,5X, 'C  =',F16.4/ 

2  '  DLJ  = ' , F16 . 8 , 5X, ' DLJK= ' , F16 . 8 , 5X, ' DLK  =',F16.8/ 

3  '  DSJ  = ' , F16 . 8 , 5X, ' DSK  =',F16.8//) 

DO  100  1=1,3 

100  LTYPE ( I ) = • FALSE 

IF  (IITY.GE. 100)  LTYPE (1)=. TRUE. 

IF  (MOD(IITY, 100) .GE. 10)  LTYPE ( 2 )=. TRUE . 

IF  (MOD ( IITY , 10 )  .GE.  1)  LTYPE ( 3 )=. TRUE . 

WRITE (6, 6030)  LTYPE , JMIN , JMAX 

6030  FORMAT (//'  SELECTION  RULES  ',3L4/ 

1  '  J  RANGE  '  ,14,  '  TO  '  ,14) 

WRITE (6,6040)  KTAU , FMIN , FMAX , FNUZ , FNU1 , FNU2 , FNU3 , 

1  STMIN, TT,NBINS, WIDTH , BS 

6040  FORMAT (//'  MAX.  DELTA  TAU  ',14/ 

1  '  FREQ.  RANGE  ' ,F10.4,'  TO  ',F10.4,'  CM-1' 

2  '  BAND  ORIGIN  ' ,F10.4,'  CM-1'/ 

7  '  NU  #1  ' , F10 . 4 , '  CM-1'/ 

8  '  NU  #2  ' , F10 . 4 , '  CM-1'/ 

9  '  NU  #3  ' , F10 . 4 , '  CM-1 '/ 

3  '  MIN.  STRENGTH  ',E12.3/ 

4  '  TEMPERATURE  ' , F7 . 2/ 

5  '  NBINS  ',15/ 

6  '  WIDTH  '  ,  F10 . 6/ 

7  '  BAND  STRENGTH  *,F10.4///) 


ooo  n  oo  o  n  n  n  n  o\  noo  no  ooo  ooo  ooo  non  no 


SET  LINE  COUNTER  TO  ZERO 
NLINE  =  0 

SET  INITIAL  TOTAL  BAND  INTENSITY  =  0 
TOTINT=0 . 0 

DEFINE  FCT  =  FACTOR  FOR  BOLTZMANN  CALCULATIONS 
FCT=-1 . 438786/TT 

CALCULATE  VIBRATIONAL  PARTITION  FUNCTION 

VIBPF=l./( (l.-DEXP(FCT*FNUl) ) * ( 1 . -DEXP ( FCT*FNU2 ) ) * 
1  ( 1 . “DEXP ( FCT*FNU3 ) ) ) 

CALCULATE  ROTATIONAL  PARTITION  FUNCTION 

ROTPF=(TT**3/ (PAR (16) *PAR ( 17 ) *PAR ( 18 ) ) ) **. 5 

DO  9000  J  =JMIN , JMAX 

PRINT  J  ON  SCREEN  TO  TRACK  PROGRESS  OF  CALCULATION 

PRINT  6010 ,J 
310  FORMAT  (IX, 15) 

******  r  BRANCH  TRANSITIONS 

IF  (J.EQ.JMIN)  GO  TO  4000 
FLJIN=1 , D0/DBLE ( J) 

JLO  =  J-l 

DO  2000  NTYPE=1 , 3 

IF  ( .NOT. LTYPE (NTYPE) )  GO  TO  2000 


DO  2100  ISYMF=1 , 4 

ISYMI  =  ISEL( ISYMF, NTYPE, 1) 

CALL  ECALC ( 0 , ISYMF , J , HD , HF , EL, NHF , KMAXF , MAXT) 

CALL  ECALC ( 1 , ISYMI , JLO , HDLO , HFLO , ELLO , NHI , KMAXI , MAXTLO ) 
IF  (NHF*NHI.EQ.O)  GO  TO  2100 

******  CALCULATE  DIRECTION  COSINE 


C 


CALL  CFACTR ( 1 , NTYPE , NHF , NHI , KMAXF , KMAXI , J , CFAC , CFACP) 


c******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT-MAXTLO 
NSMAX  =  (KTAU-NTDEL+100)/4-25 
NSMIN  =  2  5- ( KTAU+NTDEL+100 ) / 4 
C 
C 

DO  2200  ISHIFT=NSMIN, NSMAX 
IF  (ISHIFT. GE. 0)  GO  TO  2180 
N  =  MINO (NHF+ ISHIFT , NHI ) 

IF  (N.LT.l)  GO  TO  2200 
NPLUSF  =  -ISHIFT 
NPLUSI  =  0 
GO  TO  2210 
C 

2180  N  =  MINO (NHF , NHI -ISHIFT) 

IF  (N.LT.l)  GOTO  2200 
NPLUSF  =  0 
NPLUSI  =  ISHIFT 
2210  CONTINUE 
C 

C******  CALCULATE  INDIVIDUAL  TRANSITIONS 
C 

DO  2300  1=1, N 
INOF  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL (INOF) 

ELI  =  ELLO (INOI) 

FREQ  =  (ELF-ELI ) /I . +FNUZ 
C 

C  TEST  WHETHER  LINE  POSITION  IS  WITHIN  SPECIFIED  FREQUENCY 

C  RANGE.  IF  SO,  CONTINUE  WITH  INTENSITY  CALCULATION; 

C  IF  NOT,  CALCULATE  NEXT  LINE. 

C 

IF  ( FREQ. GT.FMAX. OR. FREQ. LT.FMIN)  GO  TO  2300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2* (1-INOF) 

KCF  =  J+l-KAF-ISYMF+2* (ISYMF/2) 

KAI  =  KMAXI+2* ( 1-INOI ) 

KCI  =  JLO+l-KAI-ISYMI+2* (ISYMI/2) 

C 

C******  CALCULATE  EIGEN  VECTORS 
C 

CALL  EIVEC ( NHF , HD , HF , ELF , TVF ) 

CALL  EIVEC (NHI , HDLO , HFLO , ELI , TVI ) 

C 

C******  CALCULATE  LINE  STRENGTH 
C 

FAC  =  0.0D0 
IPLUS  =  0 

IF  ( KMAXF-KMAXI . EQ . 2 )  IPLUS  =  1 


NN  =  MINO (NHF-IPLUS , NHI) 

DO  2520  11=1, NN 

2520  FAC  =  FAC+CFAC (II) *TVF ( II+IPLUS) *TVI (II) 

IF  (NTYPE.EQ.l)  GOTO  2540 

FACP  =  0.0D0 

NN  =  MINO (NHF-1, NHI) 

IF  (NN.LT.l)  GO  TO  2512 
DO  2530  11=1, NN 

2530  FACP  =  FACP+CFACP ( II ) *TVF ( II+l ) *TVI ( II ) 

2512  IF  (NTYPE.EQ.2)  FACP  =  -FACP 
FAC  =  FAC+FACP 
2540  STRE  =  FAC**2*FLJIN 
C 

C  MULTIPLY  DIRECTION  COSINES  FACTOR  BY  ROTATIONAL  BOLTZMANN 

C  FACTOR  OVER  ROTATIONAL  PARTITION  FUNCTION 

C 

STRE= (STRE*DEXP ( -ELF/ ( 0 . 695*TT) ) ) /ROTPF 
C 

C  ASSIGN  NUCLEAR  STATISTICS  FACTOR  BASED  ON  IDENTITY  OF 

C  ISOTOPE  AND  WHETHER  LEVEL  IS  ODD  OR  EVEN 

C  H20  1:3  FOR  EVEN: ODD 

C  D20  6:3  FOR  EVEN: ODD 

C  HDO  1  FOR  ALL  LEVELS 

C 

KK=IABS (KAF-KCF) 

IF ( ISOTOPE. EQ. 161)  THEN 
STNUC= ( 1 . +2 . * (REAL (MOD (KK, 2 ) ) ) ) 

ELSEIF( ISOTOPE. EQ. 2 62)  THEN 
STNUC=  ( 6  .  -3" .  *  (REAL (MOD (KK,  2 )  )  )  ) 

ELSE 

STNUC=1 . 0 
ENDIF 
C 

C  STIMULATED  EMISSION  IS  NOT  USED  IN  EMISSION  CALCULATION 

C 

C  STIMEM= ( 1 . -DEXP (FCT*FREQ) ) 

C 

C  CALCULATE  STRENGTH  WHICH  INCLUDES  ROTATIONAL  AND  NUCLEAR 

C  STATISTICS  FACTORS.  THIS  WILL  BE  USED  LATER  TO  NORMALIZE 

C  FOR  COMPARISON  OF  DIFFERENT  ISOTOPES.  TOTINT  IS  SUM  OF 

C  STRENGTH  FOR  EACH  LINE. 

C 

STRENGTH=STNUC *  S  TRE 
C 

C  CALCULATE  FINAL  INTENSITY  BY  MULTIPLYING  BY  FREQUENCY* *4 , 

C  BAND  STRENGTH,  VIBRATIONAL  BOLTZMANN  FACTOR,  AND  DIVIDING 

C  BY  VIBRATIONAL  PARTITION  FUNCTION  AND  BAND  CENTER. 

C 

STRE= (STRENGTH* (FREQ**4) *BS *DEXP ( - ( FNUZ ) / ( 0 . 69 5 *TT) )  )  ' 

1  (VIBPF*FNUZ ) 

C 

C******  STORE  THE  CALCULATED  TRANSITION  IF  STRE  EXCEEDS  THE 
C  MINIMUM  STRENGTH  CRITERION. 


c 

IF  (STRE.GE.STMIN)  CALL  STORE (NLINE , J , KAF, KCF, JLO , KAI , 

1  KCI , FREQ , STRE , NTYPE ) 

2300  CONTINUE 
2200  CONTINUE 
2100  CONTINUE 
2000  CONTINUE 
C 
C 

C******  p  BRANCH  TRANSITIONS  -  P  BRANCH  TRANSITION  FREQUENCIES 
C  AND  STRENGTHS  ARE  CALCULATED  BY  THE  SAME  METHOD  USED 

C  ABOVE  FOR  THE  R  BRANCH  TRANSITIONS. 

C 

4000  CONTINUE 

C 

C 

IF  (J.EQ.JMIN)  GO  TO  4200 
FLJIN  =  1 . DO/DBLE ( J) 

JLO  =  J-l 
C 
C 

DO  3000  NTYPE=1 , 3 

IF  ( .NOT. LTYPE( NTYPE ) )  GO  TO  3000 
C 

DO  3100  ISYMF=1,4 

ISYMI  =  ISEL(ISYMF, NTYPE, 1) 

CALL  ECALC ( 1 , IS YMF , J , HD , HF , EL , NHF , KMAXF , MAXT ) 

CALL  ECALC ( 0 , ISYMI , JLO , HDLO , HFLO , ELLO , NHI , KMAXI , MAHTLO ) 
IF  (NHF*NHI . EQ . 0 )  GO  TO  3100 
C 

C******  CALCULATE  DIRECTION  COSINES 
C 

CALL  CFACTR ( 1 , NTYPE , NHF , NHI , KMAXF , KMAXI , J , CFAC , CFACP ) 

C 

C******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT -MAXT  LO 
NSMAX  =  (KTAU-NTDEL+100 ) /4-25 
NSMIN  =  25- (KTAU+NTDEL+100 ) /4 
C 
C 

DO  3200  ISHIFT=NSMIN , NSMAX 
IF  (ISHIFT. GE. 0)  GO  TO  3180 
N  =  MINO (NHF+ ISHIFT , NHI ) 

IF  (N.LT.l)  GOTO  3200 
NPLUSF  =  -ISHIFT 
NPLUSI  =  0 
GOTO  3210 
C 

3130  N  =  MINO (NHF, NHI-ISHIFT) 

IF  (N.LT.l)  GOTO  3200 
NPLUSF  =  0 
NPLUSI  =  ISHIFT 
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3210  CONTINUE 
C 

C******  CALCULATE  INDIVIDUAL  TRANSITIONS 
C 

DO  3300  1=1, N 
INOF  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL (INOF) 

ELI  =  ELLO(INOI) 

FREQ  =  (ELI -ELF) / 1 . +FNUZ 

IF  ( FREQ. GT.FMAX. OR. FREQ. LT.FMIN)  GO  TO  3300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2* ( 1-INOF) 

KCF  =  J+l-KAF-ISYMF+2* (ISYMF/2) 

KAI  =  KMAXI-J-2*  ( 1-INOI) 

KCI  =  JLO+l-KAI-ISYMI+2* (ISYMI/2) 

C 

C******  CALCULATE  EIGEN  VECTORS 
C 

CALL  EIVEC (NHF, HD, HF, ELF, TVF) 

CALL  EIVEC (NHI , HDLO , HFLO , ELI , TVI ) 

C 

C******  CALCULATE  LINE  STRENGTH 
C 

FAC  =  O.ODO 
I PLUS  =  0 

IF  ( KMAXF-KMAXI . EQ . 2 )  IPLUS  =  1 
NN  =  MINO (NHF-IPLUS , NHI) 

DO  3520  11=1, NN 

3520  FAC  =  FAC+CFAC (II) *TVF ( II+IPLUS ) *TVI ( II ) 

IF  (NTYPE.EQ.l)  GOTO  3540 

FACP  =  O.ODO 

NN  =  MINO (NHF-i, NHI) 

IF  (NN.LT.l)  GO  TO  3512 
DO  3530  11=1, NN 

3530  FACP  =  FACP+CFACP (II) *TVF ( II+l) *TVI ( II ) 

3512  IF  (NTYPE.EQ.2)  FACP  =  -FACP 
FAC  =  FAC+FACP 
3540  STRE  =  FAC**2*FLJIN 

STRE= (STRE*DEXP ( -ELF/ ( 0 . 69  5  *TT) ) )/ROTPF 
KK=IABS (KAF-KCF) 

IF(ISOTOPE. EQ. 161)  THEN 
STNUC= ( 1 . +2 . * (REAL (MOD ( KK, 2 ) ) ) ) 

ELSEIF ( ISOTOPE. EQ. 262)  THEN 
STNUC= (6.-3.* (REAL (MOD ( KK , 2 ) ) ) ) 

ELSE 

STNUC=1 . 0 
ENDIF 

C  STIMEM= ( 1 . -DEXP ( FCT*FREQ) ) 

STRENGTH=STNUC*STRE 

STRE= ( STRENGTH* ( FREQ* *4 ) *BS*DEXP ( - ( FNUZ ) / ( 0 . 695*TT)  )  ) , 


1  ( VIBPF*FNUZ ) 

C 

C******  STORE  THE  CALCULATED  TRANSITION 
C 

IF  (STRE.GE.STMIN)  CALL  STORE (NLINE , JLO , KAI , KCI , J , KAF , 
1  KCF , FREQ , STRE , NTYPE ) 


c 

3300 

CONTINUE 

3200 

CONTINUE 

3100 

CONTINUE 

3000 

C 

C 

CONTINUE 

4100 

C 

C 

CONTINUE 

C******  Q  BRANCH  TRANSITIONS  -.Q  BRANCH  TRANSITION  FREQUENCIES 
C  AND  LINE  STRENGTHS  ARE  CALCULATED  BY  THE  SAME  METHOD 

C  USED  ABOVE  FOR  THE  P  AND  R  BRANCHES. 

C 

IF  (J.EQ.O)  GOTO  4200 

FLJIN  =  DBLE (J+J+l) /DBLE ( J* ( J+l) ) 

C 

C 

DO  7000  NTYPE=1,3 

IF  ( .NOT. LTYPE( NTYPE ) )  GO  TO  7000 
C 
C 

DO  7100  ISYMF=1,4 

ISYMI  »  ISEL (IS YMF, NTYPE, 0) 

CALL  ECALC ( 0 , IS YMF , J , HD , HF , EL , NHF , KMAXF , MAXT ) 

CALL  ECALC ( 1 , ISYMI , J , HDLO , HFLO , ELLO , NHI , KMAXI , MAXTLO ) 

IF  (NHF*NHI . EQ. 0)  GO  TO  7100 
C 

C  *****  *  CALCULATE  DIRECTION  COSINES 
C 

CALL  CFACTR ( 0 , NTYPE , NHF , NHI , KMAXF , KMAXI , J , CFAC , CFACP) 

C 

C******  RANGE  OF  ISHIFT 
C 

NTDEL  =  MAXT-MAXTLO 
NSMAX  *  (KTAU-NTDEL+100) /4-25 
NSMIN  =  25- (KTAU+NTDEL+100) /4 
C 
C 

DO  7200  ISHIFT=NSMIN, NSMAX 
C 

IF  (ISHIFT. GE. 0)  GO  TO  7320 
N  =  MINO (NHF+ ISHIFT , NHI ) 

IF  (N.LT.l)  GO  TO  7200 
NPLUSF  =  -ISHIFT 
NPLUSI  =  0 
GO  TO  7330 


7320  N  =  MINO (NHF, NHI-ISHIFT) 

IF  (N.LT.l)  GO  TO  7200 
NPLUSF  =  0 
NPLUSI  =  ISHIFT 
C 

c******  INDIVIDUAL  TRANSITIONS 
C 

7330  DO  7300  1=1, N 

INOF  =  I+NPLUSF 
INOI  =  I+NPLUSI 
ELF  =  EL (INOF) 

ELI  =  ELLO(INOI) 

FREQ  =  ( ELF-ELI ) / 1 . +FNUZ 

IF  ( FREQ. GT.FMAX. OR. FREQ. LT.FMIN)  GO  TO  7300 
C 

C******  LEVEL  NAMES 
C 

KAF  =  KMAXF+2  * ( 1-INOF) 

KCF  =  J+l-KAF-MOD ( IS YMF , 2 ) 

KAI  =  KMAXI+2* (1-INOI) 

KCI  =  J+l-KAI-MOD ( ISYMI , 2 ) 

C 

C******EIGEN  VECTORS 
C 

CALL  EIVEC (NHF , HD, HF, ELF, TVF) 

CALL  EIVEC (NHI , HDLO , HFLO , ELI , TVI ) 

C 

C******LINE  STRENGTH 
C 

FAC  =  0.0D0 

NN  =  MINO (NHI, NHF) 

DO  7520  11=1, NN 

7520  FAC  =  FAC+CFAC ( II ) *TVF (II) *TVI (II) 

IF  (NTYPE.EQ.l)  GO  TO  7540 
C 

FACP  =  0.0D0 

IF  (KMAXF.LT.J)  GO  TO  7510 
NN  =  MINO (NHF-1, NHI) 

IF  (NN.LT.l)  GO  TO  7512 
DO  7514  11=1, NN 

7514  FACP  =  FACP+CFACP ( II ) *TVF (II+l) *TVI ( II ) 

GOTO  7512 
C 

7510  NN  =  MINO (NHF, NHI-1) 

IF  (NN.LT.l)  GO  TO  7512 


DO  7516  11=1, NN 

7516  FACP  =  FACP+CFACP (II) *TVF (II) *TVI (II+l) 
C 

7512  IF  (NTYPE.EQ.3)  FACP  =  -FACP 
FAC  =  FAC+FACP 
C 

7540  STRE  =  FAC**2*FLJIN 
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STRE= ( STRE*DEXP ( -ELF/ (0.695  *TT) ) )/ROTPF 
KK=IABS (KAF-KCF) 

IF ( ISOTOPE. EQ. 161)  THEN 
STNUC= ( 1 . +2 . * (REAL (MOD ( KK, 2 ) ) ) ) 

ELSEIF ( ISOTOPE . EQ . 262 )  THEN 
STNUC= (6.-3.* (REAL (MOD ( KK , 2 ) ) ) ) 

ELSE  * 

STNUC=1. 0 
ENDIF 

C  STIMEM= ( 1 . -DEXP ( FCT*FREQ) ) 

STRENGTH=STNUC*STRE 

STRE= (STRENGTH* (FREQ* *4) *BS*DEXP ( - ( FNUZ ) / ( 0 . 695*TT) ) )/ 

1  (VIBPF*FNUZ) 

C******STORE  THE  CALCULATED  TRANSITION 
C 

IF  (STRE.GE.STMIN)  CALL  STORE (NLINE , J , KAF , KCF , J , KAI , 

1  KCI , FREQ, STRE , NTYPE) 

C 

7300  CONTINUE 
7200  CONTINUE 
7100  CONTINUE 
7000  CONTINUE 
C 
C 

4200  CONTINUE 

C 

C 

9000  CONTINUE 

C 

C 

C******  SORT  AND  PRINTOUT  THE  CALCULATED  SPECTRUM 
C 

C  PRINT  OUT  THE  TOTAL  NUMBER  OF  LINES  MEETING  THE  STRENGTH 

C  AND  FREQUENCY  RANGE  CRITERIA. 

C 

WRITE (6,6013) ' NLINE 

6013  FORMAT (1 OX, 'TOTAL  NUMBER  OF  LINES  =  ’,16/) 

C 

C  ARRANGE  THE  STORED  TRANSITIONS  IN  ORDER  OF  INCREASING 

C  FREQUENCY . 

C 

CALL  SORT (NLINE) 

C 

C  PRINT  OUT  THE  TRANSITION  TYPE  AND  QUANTUM  STATE  LABELS, 

C  THE  TRANSITION  FREQUENCY,  THE  LINE  STRENGTH  AND  THE 

C  RELATIVE  STRENGTH.  DISABLE  THIS  SUBROUTINE  CALL  TO  SUPPRES 

C  PRINTING. 

C 

C  CALL  PRINT (NLINE) 

C 

C  PRINT  OUT  THE  SUM  OF  STRENGTH  WHICH  INCLUDES  ONLY  THE 

C  DIRECTION  COSINE  AND  THE  ROTATIONAL  TERMS. 

C 


ooon  noon 


WRITE  (6,6017)  TOTINT 

6017  FORMAT (1 OX, 'TOTAL  BAND  INTENSITY  =  *,1PE10.4/) 

SORT  THE  STORED  TRANSITIONS  INTO  FREQUENCY  BINS  OF  WIDTH 
EQUAL  TO  THE  FREQUENCY  RANGE  DIVIDED  BY  NBINS . 

CALL  BINSORT (NLINE) 

APPLY  A  GAUSSIAN  LINE  SHAPE  OF  HALFWIDTH  =  "WIDTH"  TO  THE 
SUMMED  INTENSITY  IN  EACH  BIN. 

CALL  SHAPE (NLINE) 

C 

STOP 

END 


:;$l 
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INTEGER  FUNCTION  ISEL( ISYM, NTYPE , IQR) 


DIMENSION  ITBL (4,3,2) 

DATA  ITBL/  2, 1,4,3,  3, 4, 1,2,  4, 3,2,1, 

1  1,2, 3,4,  4, 3, 2,1,  3, 4, 1,2/ 

ISEL  =  ITBL(ISYM, NTYPE, IQR+1) 

RETURN 

END 

SUBROUTINE  CFACTR ( IQR , NTYPE , NHF , NHI , KMAXF , KMAXI , J , 
1CFAC , CFACP) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  CFAC(50) ,CFACP(50) 

C 

C  THIS  ROUTINE  CALCULATES  THE  DIRECTION  COSINES 

C  J  NUMBER  OF  I  STATE  SHOULD  BE  LESS  THAN  OR  EQUAL  TO 

C  THAT  OF  F  STATE 

C 

N  =  MINO (NHF , NHI ) 

C 

IF  (IQR.EQ.O)  GO  TO  1000 
C 

C******  P , R  BRANCH 

C 

KP  =  KMAXI +2 

IF  (NTYPE. NE.l)  GO  TO  910 
C 

C******  A  TYPE  TRANSITION 
C  - 

JJ  =  J*J 
DO  920  1=1, N 
K  =  KP-2  *1 

920  CFAC(I)  =  DSQRT (DBLE ( JJ-K*K) ) 

RETURN 

C 

C******  B, C  TYPE  TRANSITIONS 
C 

910  DO  940  1=1, N 
K  =  KP-2*I 
W  =  (J+K) * (J+K+l) 

IF  (K.EQ.O)  W  =  W+W 
940  CFAC(I)  =  0 . 5D0*DSQRT (W) 

C 

N  =  MINO (NHF-1, NHI) 

IF  (N.LE.O)  RETURN 
KP  =  KP- 1 
DO  950  1=1, N 
K  =  KP-2  *1 
W  =  (J-K)  *  (J-K-l) 

IF  (K.EQ.O)  W  =  W+W 


950 


CFACP(I)  =  0. 5D0*DSQRT(W) 
RETURN 


GO  TO  1160 


C******  Q  BRANCH 
C 

1000  KP  =  KMAXF+2 

IF  (NTYPE.NE. 1) 

C 

C******  A  TYPE  TRANSITION 
C 

DO  1170  1=1, N 
1170  CFAC(I)  =  KP-2*I 
RETURN 
C 

C******  B, C  TYPE  TRANSITIONS 
C 

1160  JJ  =  J*(J+1) 

IF  (KMAXF.LT. J)  GO  TO  1190 
KP  =  KP-1 

NP  =  MINO (NHF-1 , NHI ) 

GO  TO  1165 
NP  =  MINO (NHF,NHI-1) 


1190 

C 

1165 


1120 

C 


1125 


DO  1120  1=1, N 
K  =  KP— 2*1 
W  =  JJ-K* (K+l ) 

IF  (K.EQ.O)  W  =  W+W 
CFAC(I)  =  0. 5D0*DSQRT(W) 

IF  (NP.LT.l)  RETURN 
KP  =  KP-1 
DO  1125  1=1, NP 
K  =  KP-2  * I 
W  =  JJ-K* (K+l) 

IF  (K.EQ.O)  W  =  W+W 
CFACP(I)  =  0 . 5DQ*DSQRT (W) 

RETURN 

END 

SUBROUTINE  ECALC(IUL, ISYM , J , HD , HF , EL, NH , KMX , MT) 
IMPLICIT  REAL* 8 ( A-H , O-Z ) 

COMMON  /CON/  PAR (30)  /CONST/  P(15) 

DIMENSION  HD ( 50 ) ,HF(50) ,EL(50) ,HF1(50) 


C 

100 

110 

c 

c 


IF  (J.LE.98) 
NH  =  0 
RETURN 


GO  TO  100 


DO  110  1=1,15 

P(I)  =  PAR ( I+15*IUL) 


CALL  MATRIX ( ISYM, HD , HF, NH, KMX , J ) 
IF  (NH.EQ.O)  RETURN 
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EL(NH)  =  HD(NH) 

DO  200  1=1 , NH-1 
HF1 (I)  =  HF ( I ) **2 
EL{I)  =  HD ( I ) 

200  CONTINUE 
C 

IF  (NH.NE.l)  CALL  DTSM (NH , EL, HF1 ) 

MT  =  3*J-ISYM+l-(J/2) *4 
IF  (MT.GT.J)  MT  =  MT-4 
RETURN 
END 

SUBROUTINE  STORE ( NLINE , J , KAF , KCF , JL,  KAI , KCI , FREQ , 
1STRE , NTYPE) 

IMPLICIT  REAL*8 ( A-H, 0-Z) 


THIS  ROUTINE  STORES  THE  FREQUENCY,  STRENGTH,  AND  LEVEL 
NAME  OF  ALL  TRANSITIONS  MEETING  THE  STRENGTH  AND 
FREQUENCY  RANGE  CRITERIA. 


COMMON  /TBL/  NLEV ( 10000 ), FRQ ( 10000 ), ST ( 10000 ) 
COMMON/ INTEN/TOTINT, STRENGTH 


THE  LINE  COUNTER  NLINE  IS  ADVANCED  AS  EACH  ADDITIONAL 
LINE  IS  STORED. 

NLINE  =  NLINE+1 

NO  MORE  THAN  10,000  LINES  MAY  BE  STORED 

IF  (NLINE. GT. 10000)  RETURN 
FRQ (NLINE)  =  FREQ 
ST (NLINE)  =  STRE 

NLEV (NLINE) =  ( (J*100+KAF) *2+ (KCF+KAF-J) ) *100000 

1  +  ( (JL*100+KAI) *2+(KCI+KAI-JL) ) *5 

2  + NTYPE 

STRENGTH  IS  SUMMED  TO  GIVE  TOTINT. 

TOTINT=TOTINT+STRENGTH 

RETURN 

END 

SUBROUTINE  SORT (NLINE) 

THIS  ROUTINE  SORTS  THE  STORED  TRANSITIONS  BY  FREQUENCY 
AND  STORES  THEM  AGAIN  IN  ORDER  OF  INCREASING  WAVENUMBER. 

IMPLICIT  REAL*8 (A-H , O-Z ) 

COMMON  /TBL/  NLEV ( 10000 ), FRQ ( 10000 ), ST ( 10000 ) 

COMMON/ INTEN/TOTINT , STRENGTH 

LOGNB2  =  INT ( ( ALOG ( FLOAT (NLINE) ) /ALOG (2.0) ) +1 . E-5 ) 
M=NLINE 
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DO  49  NN  =  1 , LOGNB2 
M=M/2 

K  =  NLINE  -  M 
DO  44  J  =  1,K 
I  =  J 
3  CONTINUE 

L  =  I+M 

IF  (FRQ(L) .LT.FRQ(I) )  THEN 
F  =FRQ ( I ) 

S  =ST ( I ) 

N  =NLEV (I) 

FRQ ( I ) =FRQ (L) 

FRQ (L) =F 
ST(I)=ST(L) 

ST (L) =S 

NLEV ( I ) =NLEV ( L) 

NLEV(L) =N 
I  =  I-M 

IF  (I.GE.l)  GOTO  3 
ENDIF 

44  CONTINUE 

49  CONTINUE 

C 

C  THIS  LOOP  NORMALIZES  EACH  ISOTOPE  SO  THAT  THE  SUM  OF 

C  TOTINT  IS  EQUAL  TO  1.  THIS  IS  DONE  BY  DIVIDING  THE 

C  STRENGTH  OF  EACH  LINE  BY  TOTINT. 

C 

DO  8005  1=1, NLINE 
ST (I) =ST (I) /TOTINT 
8005  CONTINUE 

RETURN 
END 

SUBROUTINE  PRINT (NLINE) 

U 

C  THIS  ROUTINE  PRINTS  THE  TRANSITION  LABELS,  FREQUENCY, 

C  STRENGTH,  AND  RELATIVE  STRENGTH  FOR  EACH  LINE. 

C 

IMPLICIT  REAL* 8 ( A-H , O-Z ) 

COMMON  /TBL/  NLEV ( 10000 ), FRQ ( 10000 ), ST ( 10000 ) 
DIMENSION  INTYPE (3) , ITYP(3) 

DATA  INTYPE/ 1HA , 1HB , 1HC/ ,  ITYP/1HP, 1HQ, 1HR/ 

C 

C 

C  SMAX,  USED  TO  CALCULATE  RELATIVE  STRENGTH,  IS 

CINITIALIZED  TO  0. 

C 

SMAX  =  0.0D0 

DO  100  1=1, NLINE 

IF  (SMAX. LT.ST(I) )  SMAX=ST(I) 

100  CONTINUE 

C 

C 

DO  200  1=1, NLINE 
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N  =  NLEV ( I ) 

NTYPE=  MOD (N, 5) 

N  =  N/5 
M  =  MOD(N, 20000) 

N  =  N/20000 

J  =  N/200 

KAF  =  MOD (N/2 , 100) 

KCF  =  J-KAF+MOD (N , 2 ) 

JLO  =  M/200 

KAI  =  MOD (M/2, 100) 

KCI  =  JLO-KAI+MOD (M, 2 ) 
RST  =  ST (I) /SMAX 


IQR  =  1 
IF  (J.EQ.JLO) 
IF  (J.GT.JLO) 
C 
C 

WRITE (6, 1000) 

1 

200  CONTINUE 
RETURN 


IQR=2 

IQR=3 


INTYPE (NTYPE) ,ITYP(IQR) , J , KAF , KCF , JLO , 
KAI , KCI , FRQ ( I ) , ST ( I ) , RST 


1000  FORMAT (5X,A1, 3X,A1, 3X, 313 , 5X, 313 , F16. 5, 3X, 1PE10. 3 , 0PF8 
END 

SUBROUTINE  EIVEC (NS , A, B, ENERGY , T) 

C  CALCULATION  OF  EIGENVECTOR 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

C  A  DIAGONAL  ELMENTS  OF  TRIDIAGONAL  MATRIX  (  LARGEST 

CTO  SMALLER) 

C  B  OFF  DIAGONAL  ELEMENTS 

DIMENSION  A(50) ,B(50) ,T(50) 

IF  (NS.GT.l)  GO  TO  100 

T ( 1 )  =  l.ODO 

RETURN 

100  WORK  =DABS (A (1) -ENERGY) 

NU=1 


.4) 


DO  10  1=2, NS 
WORKP=DABS ( A ( I ) -ENERGY) 

IF (WORKP . GT . WORK)  GO  TO  10 
WORK=WORKP 


NU=I 


10  CONTINUE 

C  CALCULATE  RATIO  OF  EIGENVECTOR 

IF (NU . GT . 1 )  T(l) =-B(l)/ (A(l) -ENERGY) 

I F ( NU . LE . 2 )  GO  TO  18 
DO  17  J=2 , NU-1 

17  T(J)=-B(J)/ (B(J-l) *T(J-1) +A(J) -ENERGY) 

13  IF(NU.LT.NS)  T (NS ) =-B (NS-1 ) / ( A (NS ) -ENERGY) 

IF (NU.GE.NS-1)  GO  TO  20 
K=NS-1-NU 


DO  19  J=1 , K 
JA=NS-J 
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19  T  ( JA)  =-B  ( JA-1)  /  (  3  ( JA)  *T  ( JA-rl ) -A  ( JA)  -ENERGY) 

C  CALCULATE  EIGEN-VECTOR 

20  T ( NU ) = 1 . 0 
T2=T(NU) *  *  2 

IF (NU . EQ . 1 )  GO  TO  22 
DO  21  J  = 1 , NU - 1 
JA=NU-J 

T ( JA) =T(JA+1) *T(JA) 

21  T2=T2+T(JA) **2 

22  IF (NU . EQ . NS )  GO  TO  24 

DO  23  J=NU+1,NS 

T(J) =T ( J-l ) *T(J) 

23  T2=T2+T(J) **2 

24  T2=l . 0/DSQRT (T2 ) 

DO  25  J=1 , NS 

25  T ( J) =T ( J) *T2 
RETURN 

END 

SUBROUTINE  DTSM (NS , Q , E ) 

C  DIAGONALIZATION  OF  TRIDIAGONAL  MATRIX 

IMPLICIT  REAL* 8 ( A-H , O-Z ) 

C  Q  DIAGONAL  FROM  THE  LARGEST  ELEMENT  TO  SMALLER 

C  NS  MUST  BE  LARGER  THAN  1 

DIMENSION  Q ( 50 ) , E ( 50 ) 

C  E  (QFF-DI AGONAL  ELEMENT) **2 

PASS=0 . ID— 4 
E (NS ) =0 . 0 
C=0.0 
NR=NS 

C  SHIFT  ORIGIN 

5  RE2=DSQRT ( E ( 1 )  ) 

CM=Q ( 1 ) -RE2 

DO  6  1=2, NR 
RE1=RE2 

RE2=DSQRT ( E ( I )  ) 

CMP=Q(I) -RE1-RE2 
IF ( CM. GT . CMP) CM=CMP 

6  CONTINUE 

DO  7  1=1, NS 

7  Q(I)=Q(I)-CM 
C=C+CM 

DO  8  1=2, NR 
E(I-1)=E(I-1)/Q(I-1) 

8  Q ( I ) =Q (I) -E(I-l) 

C  REPEAT  ORTHOGONAL  TRANSITION 

NROT=0 

11  IF (NROT . GE . 10 )  GO  TO  102 

NROT=NROT+l 
Q(1)=Q(1)+E(1) 

DO  9  1=2, NR 

E  ( 1-1) =Q ( I ) *E(I-1)/Q(I-1) 

9  Q ( I) =Q (I)-E(I-1)+E(I) 

IF(E(NR-1) -PASS)  103,11,11 


5! 


;v, 


102 


12 

103 

10 

C 


20 


DO  12  1=1, NR- 1 
Q(I)=Q(I)+E(I) 

E ( I ) =E ( I ) *Q(I+1) 

GO  TO  5 
NR=NR- 1 
E (NR) =0 . 0 

IF(NR.GT.l)  GO  TO  102 

DO  10  1=1, NS 

Q(I)=Q(I)+C 

RETURN 

END 


CONSTRUCTION  OF  TRIDIAGONAL  MATRIX 
SUBROUTINE  MATRIX ( ISYM, HDISYM , HOISYM , NH , KMAX , J ) 
IMPLICIT  REAL* 8 ( A-H , O-Z ) 

COMMON/ CONST/ A ,  B ,  C ,  DU  ,  DUK ,  DLK ,  DS  J  ,  DSK ,  HU  ,  HUK ,  HLKJ  , 
1  HLK, HSJ , HSJK, HSK 

DIMENSION  HDISYM(50) ,HOISYM(50) 

KP AR= (IS  YM— 1 ) / 2 
KMAX= ( J+KPAR) /2*2-KPAR 
NH= (KMAX+2 ) /2 
IF (ISYM. EQ. 2)  NH=NH-1 
IF(NH.LE.O)  NH=0 
IF(NH.EQ.O)  RETURN 
JP=J* (J+l) 

BAVE  =  0 . 5D0* ( B+C) 

AP=A-BAVE 

BDIF  =0 . 25D0* (B-C) 

FUP  =JP 

ADD=  (  (HSJ*FUP-DSJ+HSJK)  *FUP+BDIF-DSK+HSK)  *FUP 

BJJ=(  (HU* FUP- DU)  *FUP+BAVE)  *FUP 

BJJ2=  (HUK*FUP-DUK)  *FUP+AP 

BJJ4=HLKJ*FUP-DLK 

DO  10  1=1, NH 

K=KMAX+  2-2*1 

FLKK=K*K 

HD=( (HLK*FLKK+EJJ  4 ) *FLKK+BJJ2) *FLKK+BJJ 

IF(K.NE.l)  GO  TO  10 

IF (ISYM. EQ. 3 )  HD=HD+ADD 

IF ( ISYM . EQ . 4 )  HD=HD-ADD 

HDISYM ( I) =HD 

IF(NH.LE.l)  RETURN 

BJJ=  (HSJ*FUP-DSJ+HSJK)  *FUP+BDIF-DSK+HSK 

BJJ2=HSJK*FUP-DSK+6 . 0DO*HSK 

DO  20  1=1 , NH-1 

K=KMAX+1-2*I 

FLKK=K*K 

HO=( (HSK*FLKK+BJJ2 ) *FLKK+ BJJ) *DSQRT ( DBLE ( (JP-K* ( K+ 1 ) ) 

1  * (JP-K* (K-l) ) ) ) 

IF (K. EQ. 1)  HO=HO*DSQRT ( 2 . 0D0 ) 

HOISYM ( I) =HO 
RETURN 


SUBROUTINE  SHAPE (NLINE) 


C 

C  THIS  ROUTINE  APPLIES  A  GAUSSIAN  WITH  HALFWIDTH  =  "WIDTH" 

C  TO  EACH  BIN.  AFTER  THE  LINESHAPE  APPLICATION,  IT  WRITES 

C  THE  OUTPUT  FILE  CONTAINING  THE  LIMITS  OF  THE  FREQUENCY 

C  RANGE,  NBINS ,  AND  THE  INTENSITY  FOR  EACH  BIN. 

C 

REAL  INTEN , NORM 
IMPLICIT  REAL *8 (A-H , O-Z ) 

COMMON  /TBL/  NLEV ( 10000 ), FRQ ( 10000 ), ST ( 10000 ) 

COMMON  /FREQ/  FMIN,  .'MAX 

COMMON  /BINS/  BINST ( 5000) , BINDEX ( 5000) , NBINS , WIDTH 
COMMON/OTPU77  /ACCESS  2  $ 

COMMON/  IN'T’EN/ TOTINT 

DIMENSION  INTEN (200) , SMOOTH ( 5000 ) 

CHARACTER* 2 4  ACCESS2$ 

OPEN  ( 8 , FILE=ACCESS2$ , STATUS= ’ NEW  * ,ACCESS= 

1  'SEQUENTIAL') 

WRITE  (.8,2006)  IDINT  ( FMIN)  ,  IDINT  ( FMAX) 

WRITE  (6,2006)  IDINT ( FMIN) , IDINT ( FMAX) 

WRITE (8, 2007)  NBINS 
2007  FORMAT (15) 

2006  FORMAT ( 15, 1H, , 151 

NORM=SQRT ( 2 . *  ALOG ( 2 . ) ) / ( WIDTH+ 1E-0  6 ) 

STEP=( (FMAX-FMIN) /FLOAT (NBINS) ) *NORM 
X=STEP/2 . 

2009  FORMAT (1PE10. 4) 

INTEN ( 1) =GAUSS (X) *2 . 

SUM=INTEN ( 1) 

G=INTEN(l)/2. 

M=MIN1 ( 200 . , ( 5 . 0/SNGL(STEP) )+l.) 

DO  20  J=2,M 
X=X+STEP 
INTEN (J)=-G 
G=GAUSS (X) 

INTEN (J) =G+ INTEN ( J) 

20  SUM=SUM+2 .* INTEN (J) 

BIG-0 . 0 

DO  50  1=1, NBINS 

SMOOTH ( I ) =BINST ( I ) * INTEN ( 1 ) 

DO  40  J=2 , M 
K=I+J-1 

IF (K. LE. NBINS) SMOOTH (I) =SMOOTH ( I) +BINST(K) *INTEN(J) 
K=I-J+1 

40  IF(K.GE.l)  SMOOTH ( I ) =SMOOTH ( I ) +BINST ( K) * INTEN ( J ) 

2002  FORMAT ( 14 , 5X, F10 . 4 , 5X, F10 . 4 ) 

2005  FORMAT (1PE9. 3) 

50  CONTINUE 

C 

C  TOTSMO  IS  THE  SUM  OF  THE  INTENSITES  IN  ALL  BINS  AFTER 

C  THE  GAUSSIAN  LINESHAPE  HAS  BEEN  APPLIED.  BIG  IS  THE 

C  MAXIMUM  INTENSITY  IN  ANY  OF  THE  BINS. 

C 


TOTSMO=0 . 0 
DO  49  1=1 , NBINS 
TOTSMO=TOTSMO+SMOOTH ( I ) 

IF ( SMOOTH ( I ) .GE.BIG)  BIG=SMOOTH ( I ) 

C 

C  THE  INTENSITY  IN  EACH  BIN  SMOOTH (I)  IS  WRITTEN  TO  THE 

C  OUTPUT  FILE.  * 

C 

49  WRITE ( 8 , 2005 )  SMOOTH(I) 

WRITE  (6,2011)  TOTSMO 
WRITE (6,2012)  BIG 

2012  FORMAT (1 OX, 'MAXIMUM  INTENSITY  =  ' ,1PE10.4/) 

2011  FORMAT (1 OX, 'TOTAL  BAND  INTENSITY (GAUSSIAN  SUM)  =  ', 

11PE10 . 4/ ) 

CLOSE  (8) 

RETURN 

END 

REAL  FUNCTION  GAUSS (X) 

IMPLICIT  REAL*8 ( A-H , O-Z ) 

DIMENSION  B (5) 

DATA  B/l.  3302744, -1.82 1256, 1.78 14779, -.3  56563  78,  .31938153-' 
DATA  P/. 2316419/, Q/ . 39894228/ 

T=1 ./ ( 1 . +P*X) 

PROD=0 . 0 
DO  10  1=1,5 

10  PROD=T* ( B ( I ) +PROD) 

2010  FORMAT ( 3F10 . 6) 

GAUSS=0.5-Q*PROD*EXP(-(X**2/2.0) ) 

RETURN 

END 


SUBROUTINE  BINSORT (NLINE) 


THIS  ROUTINE  SORTS  THE  STORED  TRANSITIONS  INTO  BINS  OF 
EQUAL  WIDTH  AND  SUMS  THE  INTENSITY  IN  EACH  BIN. 

IMPLICIT  REAL*8 (A-H, O-Z) 

COMMON  /TBL/  NLEV ( 10000 ), FRQ ( 10000 ), ST ( 10000 ) 

COMMON  /FREQ/  FMIN , FMAX 

COMMON  /BINS/  BINST ( 5000 ), BINDEX (5000) , NBINS , WIDTH 

THE  BIN  WIDTH  DEL  IS  CALCULATED  BY  DIVIDING  THE 
FREQUENCY  RANGE  INTO  NBINS  EQUAL  INCREMENTS. 

DEL= ( FMAX-FMIN) /FLOAT ( NBINS ) 

K=1 

DO  100  1=1, NBINS 
BINST (I) =0. 

BINDEX ( I) =  FMIN+REAL ( I*DEL) 

DO  200  J=K, NLINE 
IF (FRQ ( J) . LE . BINDEX ( I ) )  THEN 
BINST ( I )  =  BINST ( I )  +  ST(J) 

IF (J.EQ. NLINE)  GOTO  110 
ELSE 
K=J 

GO  TO  150 
ENDIF 

IF(FRQ(J) .GT. BINDEX(I) )  K=J 

CONTINUE 

CONTINUE 

CONTINUE 

DO  50  L=I+1, NBINS 
BINST (L)=0. 

CONTINUE 

RETURN 


SAMPLE  PARAMETER  FILE  FOR  C02 


ID [ 13 ] 

BU [E15 . 10 ]  DU [E15 . 10 ]  HU[E15.10] 

BUQ [E15 . 10 ]  DUQ[E15.10]  HUQ [E15 . 10 ] 

BL[E15 . 10]  DL[E15. 10]  HL[E15.10] 

JMIN [ 14 ]  JMAX [14]  FMIN [ F10 . 4 ]  FMAX[F10.4] 

FNU1 [F16 . 8 ]  FNU2 [ F16 . 8 ]  FNU3[F16.8] 

FNUZ [ F16 . 8 ]  FLOW [ F16 . 8 ]  FHIGH[F16.8] 

SV [ F16 . 8 ]  TT[F16 . 8 ] 

NBINS [ 15 ]  WIDTH [ F10 . 6 ] 


ID  =  identification  code 

BU , DU , HU  =  upper  state  constants 

BUQ , DUQ , HUQ  =  upper  state  constants  for  other  doublet  level 

BL, DL,  HL  =  lower  state  constants 

JMIN, JMAX  ■  minimum  and  maximum  values  of  J 

FMIN , FMAX  =  minimum  and  maximum  values  of  frequency 

FNUi ,  FNU2 FNU3  =  fundamental  frequencies 

FNUZ  =  band  center 

FLOW  =  lower  state  vibrational  frequency 

FHIGH  =  upper  state  vibrational  frequency 

SV  =  band  strength 

TT  =  temperature 

NBINS  =  number  of  bins 

WIDTH  =  Gaussian  halfwidth 


SAMPLE  INPUT 

FILE  C02NU3 

11 

0.38714069 

1. 32873E-7 

0 . 077E-13 

0.00 

0.00 

0.00 

0.39021817 

1. 33204E-7 

0. 055E-13 

0  20  2000. 

2500. 

1388 . 1847 

667 . 3801 

2349 . 1433 

2349.1433 

0000. 

2349.1433 

955900. 

296. 

OUTPUT  FROM  PROCRAM  C02A8S 


5/28/1987 
13:  6 

INPUT  FILE  *  C02NU3 
OUTPUT  FILE  -  C02ANU3 
ID  1  11 

UPPER  STATE  CONSTANTS 

8  =  0.38714069  D  =  1.3287E-07  H  =>  7.7000E-15 

UPPER  STATE  CONSTANTSCOTHER  DOUBLET  LEVEL) 

3  =  0.00000000  D  *  0.0000E-01  H  *  O.OOOOE-OI 

LOWER  STATE  CONSTANTS 

B  *  0.39021817  D  *  1.3320E-07  H  =  5.5000E-15 


J  RANGE  0  TO  20 

FREQUENCY  RANGE  2000.0000  TO  2500.0000 

SAND  ORIGIN  2349.1433  CM- 1 
NU1  1388.1847  CM- 1 

NU2  667.3801  CM- 1 

NU3  2349.1433  CM-1 

LOWER  STATE  VIBRATIONAL  ENERGY  *  0.0000  CM-1 

UPPER  VIBRATIONAL  STATE  ENERGY  »  2349.1433  CM-1 

SAND  STRENGTH  955900.00000000 
TEMPERATURE  296.0000 

NS  I  NS  50 00 

WIDTH  2.000000 


p 

20 

2332.3694 

1.128E+04 

8.776261 

p 

18 

2334.1569 

1 . 178E+04 

9.164751 

p 

16 

2335.9199 

1.197E+04 

9.309935 

p 

14 

2337.6586 

1.179E+04 

9.169424 

p 

12 

2339.3728 

1.120E+04 

8.713461 

p 

10 

2341.0625 

1.019E+04 

7.928887 

p 

8 

2342.7277 

8.770E*03 

6.821999 

p 

6 

2344.3685 

6.968E+03 

5.419867 

p 

4 

2345.9847 

4.846E+03 

3.769805 

p 

2 

2347.5763 

2.490E*03 

1.936939 

R 

0 

2349.9176 

1 .261E*03 

0.980531 

R 

2 

2351.4477 

3.741E+03 

2.910200 

R 

4 

2352.9531 

6.076E*03 

4.726255 

R 

6 

2354.4338 

8.164E+Q3 

6.350329 

R 

8 

2355.8899 

9.9226*03 

7.717874 

R 

10 

2357.3212 

1 .  129E+04 

8.782356 

R 

12 

2358.7277 

1.224E+04 

9.517692 

R 

14 

2360.1095 

1 .275E+04 

9.918749 

R 

16 

2361.4664 

1 .286E*04 

10.000000 

R 

18 

2362.7985 

1.259E+04 

9.792624 

2000, 

2500 

TOTAL 

BAND  I N TENS I TY( GAUSS I AN  SUM) 

=  1.8218E*05 

MAXIMUM  INTENSITY  =■  8.9935E*02 
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c 

C  PROGRAM  C02ABS . FOR 

C 

Q 

IMPLICIT  REAL*8 ( A-H , 0-2 ) 

DIMENSION  EL ( 500 ) ,EU(500) ,EROTP(500) ,EROTQ(500) , 

1  EROTR (500) , EP (500) ,EQ(500) ,ER(500) ,SP(500) , 

1  SQ(500) ,SR(500) ,EUQ(500) 

COMMON/TBL/FRQ ( 1000 ) ,ST(1000) ,NLEV(1000) 

COMMON/ FREQ/ FMI N, FMAX 

COMMON/ BINS/BINST ( 5000 ) , BINDEX ( 5000 ) ,NBINS,WIDTH 
COMMON/OTPUT/ ACCESS 2$ 

CHARACTER* 8  FILENM1$ , FILENM2$ 

CHARACTER* 2 4  ACCESS1$ , ACCESS2$ 

INTEGER* 2  I YEAR , IMONTH , I DAY , IHOUR , IMINUTE 
CALL  GETTIM( IHOUR, IMINUTE) 

CALL  GETDAT(IYEAR, IMONTH, I DAY) 

WRITE (6,*)  '  ENTER  INPUT  FILENAME?' 

READ  19 , FILENM1$ 

19  FORMAT (A8) 

WRITE (6,*)  '  ENTER  OUTPUT  FILENAME?' 

READ  19,  FILENM2  $ 

ACCESS1$='C:\ASYMRTR\'//FILENM1$// ' .DAT' 
ACCESS2$=' A: ' //FILENM2$// ' . DAT' 

OPEN ( 2 , FILE=ACCESS 1 $ , ACCESS= ’ SEQUENTIAL ' ) 

READ (2,23)  ID 
23  FORMAT (13) 

READ (2,29)  BU,DU,HU 
29  FORMAT (3E15. 10) 

READ (2, 29)  BUQ, DUQ, HUQ 
READ (2 ,29)  BL, DL, HL 
READ (2, 39)  JMIN , JMAX , FMIN , FMAX 
39  FORMAT ( 214, 2F10. 4) 

READ (2, 29)  FNU1 , FNU2 , FNU3 
READ (2, 49)  FNUZ , FLOW , FHIGH 
49  FORMAT ( 3F16 . 8 ) 

READ (2, 59)  SV,TT 
59  FORMAT (2F16 . 8 ) 

READ (2, 69)  NBINS, WIDTH 
69  FORMAT ( 15 , F10 . 6 ) 

CLOSE (2) 

WRITE (6, 1019) 

1019  FORMAT ( '  OUTPUT  FROM  PROGRAM  C02ABS ' // ) 

WRITE (6, 1025)  IMONTH , IDAY , IYEAR 
WRITE (6, 1026)  IHOUR , IMINUTE 

1025  FORMAT (IX, 12,  ' /  ’  , 12 ,  '  /  '  ,  1 4  ) 

1026  FORMAT( IX, 12 , ' : ' , 12/) 

WRITE(6, 1029)  FILENM1$ 

1029  FORMAT ( '  INPUT  FILE  =  ',A8/) 

WRITE(6, 1039)  FILENM2 $ 

1039  FORMAT ('  OUTPUT  FILE  =  ',A8/) 

WRITE(6, 1024)  ID 
1024  FORMAT (  '  ID  =  ' , 13) 


hO 


Li. 


WRITE (6, 1049)  BU,DU,HU 
1049  FORMAT ('  UPPER  STATE  CONSTANTS’/ 

1  '  B  =  ' , F16 . 8 , 3X, 1 D  =  ' , 1PE10.4 , 3X, 'H  =  ’,1PE10.4/) 

WRITE (6, 1054)  BUQ , DUQ , HUQ 

1054  FORMAT ( '  UPPER  STATE  CONSTANTS (OTHER  DOUBLET  LEVEL)  '/ 

1  '  B=  ' , F16 . 8 , 3X , ' D  =  ' , 1PE10 . 4 , 3X , ' H  =  ’,1PE10.4/) 

WRITE(6, 1059)  BL, DL, HL 
1059  FORMAT ('  LOWER  STATE  CONSTANTS'/ 

1  '  B  =  1 , F16 . 8 , 3X, ' D  =  ' , 1PE10.4, 3X, 'H  =  ’,1PE10.4/) 

WRITE (6, 1069)  JMIN , JMAX , FMIN , FMAX 
1069  FORMAT (//'  J  RANGE  ',14,'  TO  ',14/ 

1  '  FREQUENCY  RANGE  ' ,F10.4,'  TO  • ,F10.4/) 

WRITE (6, 1079)  FNUZ , FNUl , FNU2 , FNU3 


1079 

FORMAT (’  BAND  ORIGIN 

' , F10.4, ' 

CM— 1  '  / 

1 

'  NU1 

' , F10 . 4 , ' 

CM-1 '  / 

2 

'  NU2 

' , F10 . 4 , ' 

CM-1 '/ 

3 

'  NU3 

' , F10 . 4 , ' 

CM-1 '/) 

WRITE (6, 1084)  FLOW 

1084  FORMAT ('  LOWER  STATE  VIBRATIONAL  ENERGY  =  ',F10.4, 

1'  CM- 1 1 / ) 

WRITE(6, 1086)  FHIGH 

1086  FORMAT ('  UPPER  VIBRATIONAL  STATE  ENERGY  =  ',F10.4, 

1’  CM-l'/J 

WRITE (6, 1089)  SV, TT , NBINS , WIDTH 
1089  FORMAT ( '  BAND  STRENGTH  ',F16.8/ 

1  '  TEMPERATURE  ',F10.4/ 

2  '  NBINS  ' , 15/ 

3  '  WIDTH  1 , F10 . 6/// ) 

C2«l. 438786 

NLINE=0 

DO  2000  I=JMIN , JMAX 

EL(I)=BL*I*(I+1)-DL*(I**2) *( (1+1) **2 ) +HL* ( 1**3 ) * 

1( (1+1) **3 ) 

2000  CONTINUE 

DO  2500  I=JMIN , JMAX 

EU(I) =BU*I* (1+1) -DU* (1**2) * ( (1+1) **2) +HU* (1**3 ) * 
1((I=1)**3) 

2500  CONTINUE 

DO  2600  I=JMIN , JMAX 

EUQ(I) =BUQ*I* (1+1) -DUQ* (1**2) * ( (1+1) **2) + 

1  HUQ* (1**3) *( (1+1) **3) 

2600  CONTINUE 

DO  3000  I=JMIN+2 , JMAX, 2 
EROTP ( I ) =EU ( I - 1 ) -EL ( I ) 

3000  CONTINUE 

DO  3500  I=JMIN, JMAX-2 , 2 
EROTR(I) =EU(I+1) -EL (I) 

3500  CONTINUE 

VIBPFT0=1 . / ( (l.-DEXP( -C2  *FNU1/ 29  6 . )  ) * ( ( 1. -DEXP 
1  (-C2*FNU2/TT) ) **2) * ( 1 . -DEXP(-C2*FNU3/296 . ) ) ) 

VIBLZMTO=DEXP(-C2*FLOW/296. ) 

ROTPFTO=296./BL 
VIBLZM=DEXP ( -C2*FLOW/TT) 


ROTPF=TT/BL 

VIBPF=l./( ( 1 • -DEXP (-C2*FNU1/TT) )*( (l.-DEXP 
1  (-C2*FNU2/TT) ) **2 ) * ( 1 . -DEXP ( -C2  *FNU3/TT) ) ) 

S=»  (VIBLZM*SV*VIBPFTO)  /  ( VIBLZMTO*VIBPF*ROTPF) 

IF ( ID. LT. 20)  GO  TO  3900 
DO  3600  I=JMIN+2 , JMAX-2 , 2 
EROTQ ( I ) =EUQ ( I ) -EL ( I ) 

3600  CONTINUE 

DO  3700  I=JMIN+2, JMAX-2, 2 
FREQ=EROTQ ( I )  +  FNUZ 
ROTBLZM=DEXP ( (-C2*EL(I)/TT) ) 

STIMEM  =  (l.-DEXP( ( -C2*FREQ) /TT) ) 

IF (ID.EQ. 22)  THEN 
HLFAC= ( 2 . *REAL ( I ) + 1 . )  /  4  . 

ELSE 

HLFAC=(2.*REAL(I)+l.)/( (REAL(I) ) * (REAL ( I) +1 . ) ) 
ENDIF 

SQ ( I ) =S*FREQ*HLFAC*STIMEM*ROTBLZM/ FNUZ 
STRE=SQ ( I ) 

NTYPE=2000+I 

IF ( FREQ . GE . FMIN . AND . FREQ . LE . FMAX)  THEN 
CALL  STORE (NLINE , FREQ , STRE , NTYPE) 

ENDIF 

3700  CONTINUE 
3900  SMAX=0. 000001 

DO  4000  I=JMIN+2 , JMAX , 2 
FREQ=EROTP ( I ) +FNUZ 
ROTBLZM=DEXP ( ( -C2  *EL (I) ) /TT) 

STIMEM= ( 1 . -DEXP ( (-C2*FREQ) /TT) ) 

IF (ID.EQ . 11)  THEN 
HLFAC-REAL ( I ) 

ELSEIF( ID.EQ. 22)  THEN 

HLFAC=REAL ( I ) * (REAL ( I ) -1) / ( 4 . *REAL (I) ) 

ELSE 

HLFAC= (REAL( I ) +1 . ) * (REAL( I ) -1 . ) /REAL ( I ) 

ENDIF 

SP ( I ) =S*FREQ*HLFAC*STIMEM*ROTBLZM/FNUZ 
STRE=SP(I) 

NTYPE=1000+I 

IF ( FREQ . GE . FMIN . AND . FREQ . LE . FMAX )  THEN 
CALL  STORE (NLINE , FREQ , STRE , NTYPE) 

ENDIF 

4000  CONTINUE 

DO  5000  I=JMIN, JMAX-2 , 2 
FREQ=EROTR ( I ) +FNUZ 
ROTBLZM=DEXP ( (-C2*EL(I)/TT) ) 

STIMEM= ( 1 . -DEXP ( (-C2*FREQ) /TT) ) 

IF (ID. EQ. 11)  THEN 
HLFAC=REAL(I) +1. 

ELSEIF ( ID . EQ . 2  2 )  THEN 
HLFAC= ( REAL ( I ) +  2 . )/4. 

ELSE 

HLFAC= (REAL ( I ) +2 . ) *REAL ( I ) / (REAL ( I ) + 1 . ) 


.s' 

s’ 


ENDIF 

SR ( I ) =S*FREQ*HLFAC*STIMEM*ROTBLZM/FNUZ 
STRE=SR ( I ) 

NTYPE=3  000+1 

IF ( FREQ . GE . FMIN . AND . FREQ . LE . FMAX )  THEN 
CALL  STORE (NLINE, FREQ, STRE,NTYPE) 

ENDIF 

5000  CONTINUE 

CALL  SORT (NLINE) 

C  CALL  PRINT (NLINE) 

CALL  BINSORT (NLINE) 

CALL  SHAPE (NLINE) 

STOP 

END 


C  PROGRAM  C02EMS . FOR 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  EL ( 500 ) ,EU(500) ,EROTP(500) ,EROTQ(500) , 

1  EROTR ( 500 ) , EP ( 500 ) ,EQ(500) ,ER(500) ,SP(500) , 

2  SQ(500) ,SR(500) ,EUQ(500) 

COMMON/TBL/FRQ ( 1000 ) , ST (1000) ,NLEV(1000) 

COMMON/ FREQ/ FMIN , FMAX 

COMMON/BINS/BINST(5000) , BINDEX(5000) , NBINS , WIDTH 
COMMON/ OTPUT/ ACCESS2  $ 

CHARACTER* 8  FILENM1$ , FILENM2  $ 

CHARACTER* 2 4  ACCESS1$ , ACCESS2 $ 

INTEGER*2  IYEAR, IMONTH, IDAY , IHOUR, IMINUTE 
CALL  GETTIM ( IHOUR , IMINUTE ) 

CALL  GETDAT ( IYEAR , IMONTH , IDAY ) 

WRITE (6 , *)  '  ENTER  INPUT  FILENAME?' 

READ  19 , FILENM1$ 

19  FORMAT (A8) 

WRITE (6,*)  '  ENTER  OUTPUT  FILENAME?' 

READ  19,  FILENM2$ 

ACCESS1$='C:\ASYMRTR\'//FILENM1$// ' .DAT' 

ACCESS2$='A: '//FILENM2$// ' .DAT' 

OPEN ( 2 , FILE=ACCESS1$ , ACCESS= 1  SEQUENTIAL ' ) 

READ (2, 23)  ID 
23  FORMAT (13) 

READ (2,29)  BU,DU,HU 
29  FORMAT (3E15. 10) 

READ (2, 29)  BUQ, DUQ,HUQ 
READ (2,29)  BL, DL, HL 
READ (2, 3 9)  JMIN,JMAX, FMIN, FMAX 
39  FORMAT ( 214, 2F10. 4) 

READ (2, 29)  FNU1 , FNU2 , FNU3 
READ (2, 49)  FNUZ , FLOW, FHIGH 
49  FORMAT (3F16. 8) 

READ (2,59)  SV,TT 
59  FORMAT (2F16. 8) 

READ (2, 69)  NBINS, WIDTH 
69  FORMAT (15 , F10 . 6 ) 

CLOSE (2) 

WRITE (6, 1019) 

1019  FORMAT ('  OUTPUT  FROM  PROGRAM  C02EMS '//) 

WRITE (6, 1025)  IMONTH, IDAY , IYEAR 
WRITE (6, 1026)  IHOUR, IMINUTE 

1025  F0RMAT(1X, 12 ,  '/ '  , 12,  '/ '  ,14) 

1026  FORMAT(lX, 12,  ' :  '  ,12/) 

WRITE (6, 1029)  FILENM1$ 

1029  FORMAT ( '  INPUT  FILE  =  ',A8/) 

WRITE(6, 1039)  FILENM2  $ 

1039  FORMAT ( '  OUTPUT  FILE  =  ',A3/) 

WRITE(6, 1024)  ID 
1024  FORMAT (  '  ID  =  ',13) 

WRITE(6, 1049)  BU,DU,HU 
1049  FORMAT ('  UPPER  STATE  CONSTANTS'/ 

1  '  B  =  ' , F16 . 8 , 3X , ' D  =  ' , 1PE10 . 4 , 3X , ' H  =  ' , 1PE10 . 4. ) 
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WRITE (6, 1054)  BUQ,DUQ,HUQ 

1054  FORMAT ('  UPPER  STATE  CONSTANTS ( OTHER  DOUBLET  LEVEL)  ’/ 

1  '  B  =  ' , F16 . 8 , 3X , ’ D  =  ’ , 1PE10. 4 , 3X, ’H  =  ’ ,1PE10.4/) 

WRITE(6, 1059)  BL, DL, HL 
1059  FORMAT ( '  LOWER  STATE  CONSTANTS'/ 

1  '  B  =  ' , F16 . 8 , 3X , ' D  =  ' , 1PE10 . 4 , 3X, 1 H  =  ',1PE10.4/) 

WRITE (6, 1069)  JMIN , JMAX , FMIN , FMAX 
1069  FORMAT (//'  J  RANGE  ',14,'  TO  ',14/ 

1  '  FREQUENCY  RANGE  ',F10.4,'  TO  ',F10.4/) 

WRITE (6, 1079)  FNUZ , FNU1 , FNU2 , FNU3 
1079  FORMAT ( '  BAND  ORIGIN  ',F10.4,'  CM-1'/ 

1  '  NU1  ' , F10 . 4 , '  CM-1'/ 

2  '  NU2  ' , F10 . 4 , '  CM-1'/ 

3  '  NU3  ' , F10 • 4 , '  CM-1'/) 

WRITE (6, 1084)  FLOW 

1084  FORMAT ('  LOWER  VIBRATIONAL  STATE  ENERGY  =  ',F10.4, 

1'  CM-1'/) 

WRITE(6, 1086)  FHIGH 

1086  FORMAT ('  UPPER  VIBRATIONAL  STATE  ENERGY  =  ',F10.4, 

1'  CM-1'/) 

WRITE(6, 1089)  SV, TT , NBINS , WIDTH 
1089  FORMAT ('  BAND  STRENGTH  ',F16.8/ 

1  '  TEMPERATURE  ',F10.4/ 

2  '  NBINS  ',15/ 

3  '  WIDTH  ' , F10 . 6/// ) 

C2=l. 438786 

NLINE=0 

DO  2000  I=JMIN , JMAX 

EL(I)=BL*I*(I+1)-DL*(I**2) *( (1+1) **2 ) +HL* ( 1**3 ) * 

1 ( (1  +  1) **3 ) 

2000  CONTINUE 

DO  2500  I=JMIN , JMAX 

EU(I)=BU*I* (1+1) -DU* (1**2) * ( (1+1) **2) +HU* (1**3) * 

1( (1+1) **3 ) 

2500  CONTINUE 

DO  2600  I=JMIN , JMAX 

EUQ ( I ) =BUQ*I* ( 1+1) -DUQ* (1**2) * ( (1+1) **2) + 

1  HUQ* (1**3) * ( (1+1) **3) 

2600  CONTINUE 

DO  3000  I=JMIN+2 , JMAX, 2 
EROTP ( I ) =EU ( I - 1 ) -EL ( I ) 

3000  CONTINUE 

DO  3500  I=JMIN, JMAX-2 , 2 
EROTR(I)=EU(I+l) -EL ( I ) 

3500  CONTINUE 

VIBPFT0  =  1 ./ ( ( 1 • ~DEXP(-C2*  FNU1/ 2  9  6. )  ) * (  ( 1 . -DEXP 
1  ( -C2  *FNU2/2  9  6 . )  ) **2) * ( 1 . -DEXP ( -C2 *FNU3/2 9 6 . )  )  ) 
VIBLZMTO=DEXP ( -C2  *  FLOW/ 29  6 . ) 
rqtPFT0=296./BL 
VIBLZM=DEXP(-C2*FHIGH/TT) 

ROTPF=TT/BL 

VIBPF=1 . / ( ( 1 . -DEXP ( -C2  *FNU1/TT)  ) * (  ( 1 . -DEXP 
1  ( -C2  *FNU2/TT)  ) **2) * ( 1 . -DEXP ( -C2 * FNU3/TT)  )  ) 


S=(VIBLZM*SV*VIBPFTO) / (VIBLZMTO*VIBPF*ROTPF) 

IF ( ID. LT. 20)  GO  TO  3900 
DO  3600  I=JMIN+2 , JMAX-2 , 2 
EROTQ ( I ) =EUQ ( I ) -EL ( I ) 

3600  CONTINUE 

DO  3700  I=JMIN+2, JMAX-2, 2 
FREQ=EROTQ ( I )  +  FNUZ 
ROTBLZM=DEXP ( (-C 2*EL(I)/TT) ) 

C  STIMEM  =  ( 1 . -DEXP ( ( -C2  *FREQ) /TT)  ) 

IF ( ID. EQ. 22 )  THEN 
HLFAC= ( 2 . *REAL ( I ) +1 . )/4. 

ELSE 

HLFAC= ( 2 . *REAL ( I) +1 . ) / ( (REAL(I) ) * (REAL ( I ) +1 . ) ) 
ENDIF 

SQ (I) =S*FREQ**4*HLFAC*RQTBLZM/FNUZ 
STRE=SQ ( I) 

ntype=2ooo+j: 

CALL  STORE (NLINE , FREQ , STRE , NTYPE) 

3700  CONTINUE 
3900  SMAX=0. 000001 

DO  4000  I=JMIN+2 , JMAX , 2 
FREQ=EROTP ( I ) +FNUZ 
ROTBLZM=DEXP( (-C2*EL(I) )/TT) 

C  STIMEM= ( 1 . -DEXP ( ( -C2*FREQ) /TT)  ) 

IF ( ID . EQ . 11 )  THEN 
HLFAC=REAL ( I ) 

ELSE IF (ID. EQ . 22 )  THEN 

HLFAC=REAL ( I ) * ( REAL (I) -1) / (4 . *REAL ( I ) ) 

ELSE 

HLFAC= ( REAL (I)+l.)*( REAL ( I ) -1 . )/REAL(I) 

ENDIF 

SP ( I) =S*FREQ**4*HLFAC*ROTBLZM/FNUZ 
STRE=SP ( I ) 

NTYPE=1000+I 

CALL  STORE ( NLINE , FREQ , STRE , NT Y  PE ) 

4000  CONTINUE 

DO  5000  I=JMIN, JMAX-2 , 2 
FREQ=EROTR ( I ) +FNUZ 
ROTBLZM=DEXP( ( -C2 *EL ( I ) /TT) ) 

C  STIMEM= ( 1 . -DEXP ( ( -C2  *FREQ) /TT) ) 

IF(ID.EQ.ll)  THEN 
HLFAC=REAL ( I ) + 1 . 

ELSEIF ( ID . EQ . 22 )  THEN 
HLFAC= ( REAL ( I ) +  2 . ) / 4 . 

ELSE 

HLFAC= (REAL( I) +2 . ) *REAL( I ) / (REAL ( I ) +1 . ) 

ENDIF 

SR ( I ) =S*FREQ**4  *HLFAC*ROTBLZM/FNUZ 
STRE=SR ( I ) 

NTYPE  =  3  000  +  1 

CALL  STORE (NLINE , FREQ , STRE , NTYPE ) 

5000  CONTINUE 

CALL  SORT (NLINE) 
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CALL  PRINT (NLINE) 
CALL  BINSORT (NLINE) 
CALL  SHAPE (NLINE) 
STOP 
END 


C  PROGRAM  HCLEMS 

IMPLICIT  REAL*8 (A-H, O-Z) 

DIMENSION  EL ( 1000 ) ,EU(1000) ,EROTP(1000) ,EROTR(1000) 
1  EP(1000) , ER ( 1000 ) ,SP(1000) ,SR(1000) 

COMMON/TBL/FRQ( 10000) , ST (10000) 

COMMON/ FREQ/ FMIN , FMAX 

COMMON/ BINS/ BINST ( 5000 ) , BINDEX ( 5000 ) , NBINS, WIDTH 
COMMON/OTPUT/ ACCESS2$ 

COMMON/ INTEN/TOTINT, STRENGTH 
CHARACTER* 8  FILENM1$ , FILENM2$ 

CHARACTER* 2 4  ACCESS 1$ , ACCESS 2$ 

WRITE (6,*)  '  ENTER  INPUT  FILENAME?' 

READ  19 , FILENM1$ 

19  FORMAT (A8) 

WRITE (6,*)  '  ENTER  OUTPUT  FILENAME?' 

READ  19,  FILENM2  $ 

ACCESS1$='C:\ASYMRTR\'//FILENM1$// ' .DAT' 
ACCESS2$='A: ' //FILENM2$// ' .DAT' 

OPEN (2 , FILE=ACCESS1$ , ACCESS= ' SEQUENTIAL' ) 

READ (2,29)  BU,DU,HU 
29  FORMAT (3F16. 10) 

READ (2, 29)  BL,DL,HL 
READ (2, 39)  JMIN,JMAX, FMIN, FMAX 
39  FORMAT (2I4,2F10.4) 

READ (2, 49)  FNUZ 
49  FORMAT (F16. 8) 

READ (2, 59)  SV,TT 
59  FORMAT (2F16. 8) 

READ (2, 69)  NBINS, WIDTH 
69  FORMAT (15 , F10 . 6) 

CLOSE (2) 

WRITE (6, 1019) 

1019  FORMAT ( '  OUTPUT  FROM  PROGRAM  HCL ' // ) 

WRITE (6, 1029)  FILENM1$ 

1029  FORMAT (  '  INPUT  FILE  =  ' ,A8/) 

WRITE (6, 1039)  FILENM2$ 

1039  FORMAT ( '  OUTPUT  FILE  =  ’,A8/) 

WRITE(6, 1049)  BU,DU,HU 
1049  FORMAT ( '  UPPER  STATE  CONSTANTS'/ 

1  '  B  =  ' , F16 . 8 , ' D  =  ' , F16 . 8 , 1 H  = 

WRITE(6, 1059)  BL, DL, HL 
1059  FORMAT ( '  LOWER  STATE  CONSTANTS ’/ 

1  '  B  =  ' , F16 . 8 , ' D  =  ' , F16 . 8 , ' H  = 

WRITE (6, 1069)  JMIN , JMAX , FMIN , FMAX 
1069  FORMAT (//'  J  RANGE  ',14,'  TO  ',14/ 

1  '  FREQUENCY  RANGE  ',F10.4,'  TO 

WRITE (6, 1079)  FNUZ 

1079  FORMAT ( '  BAND  ORIGIN  ',F10.4,'  CM-1'/) 

WRITE(6, 1089)  SV, TT , NBINS , WIDTH 
1089  FORMAT ('  BAND  STRENGTH  ',F16.8/ 

1  '  TEMPERATURE  ',F10.4/ 

2  '  NBINS  ',15/ 

3  '  WIDTH  ' , F10 . 6/ // ) 


' , F16 . 8/ ) 
’ , F16 . 8/ ) 
'  ,F10.4/) 
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C2=l. 438786 

NLINE=0 

TOTINT=0 . 0 

DO  2000  I=JMIN , JMAX 

EL(I)=BL*I* (1+1) -DL* ( 1**2 ) *( (1  +  1) **2 ) +HL* ( 1**3  )  * 
1 ( (1  +  1) **3 ) 

1100  FORMAT ('  J=  ',14,'  ER=  • ,F10.4) 

2000  CONTINUE 

DO  2500  I=JMIN , JMAX 

EU(I) =BU*I* (1+1) -DU* (1**2) * ( (1+1) **2) +HU* (1**3) * 
1( (1+1) **3 ) 

2500  CONTINUE 

DO  3000  I=JMIN+1, JMAX 
EROTP ( I ) =EU ( I - 1 ) -EL ( I ) 

3000  CONTINUE 

DO  3500  I=JMIN , JMAX 
EROTR(I) =EU (1+1) -EL(I) 

3500  CONTINUE 

VIBPFT0=1 . / ( 1 . -DEXP ( -C2*FNUZ/296 . ) ) 
VIBL2MT0=DEXP(-C2*FNUZ/296. ) 
rOTPFT0=296./BL 
VIBLZM=DEXP ( -C2*FNUZ/TT) 

ROTPF=TT/BL 

VIBPF=1 . / ( 1 . -DEXP ( -C2*FNUZ/TT) ) 

SMAX=0 .0001 
DO  4000  I=JMIN+1, JMAX 
FREQ=EROTP ( I ) +  FNUZ 
ROTBLZM=DEXP ( ( -C2  *EU ( I ) )/TT) 

HLFAC=I 
MM=— I 

FF=l.+(-2,5599E-2) *MM+ ( 3 . 203E-4  ) *MM**2 

STRENGTH=HLFAC*ROTBLZM/ROTPF 

SP ( I ) =FREQ**4*VIBLZM*VIBPFT0*FF*STRENGTH/ 

1  (FNUZ*VIBPF) 

IF(SP(I) .GT.SMAX)  SMAX=SP ( I ) 

EP ( I ) =FREQ 
STRE=SP ( I ) 

1099  FORMAT ( 1  P  ' , 3X , 14 , 5X , F10 . 4 , 5X , 1PE9 . 3 , 5X , F10 . 6 ) 
IF ( FREQ . GE . FMIN . AND . FREQ . LE . FMAX)  THEN 
CALL  STORE (NLINE, FREQ, STRE) 

ENDIF 

4000  CONTINUE 

DO  5000  I=JMIN, JMAX-1 
FREQ=EROTR ( I ) +FNUZ 
ROTBLZM=DEXP ( ( -C2*EU ( I ) /TT) ) 

HLFAC=I+1. 

MM=I+1. 

FF=l.+(-2.5599E-2) *MM+ (3 . 203E-4 ) *MM**2 

STRENGTH=HLFAC*ROTBLZM/ROTPF 

SR ( I ) =FREQ**4  *VIBLZM*VIBPFTO*FF*STRENGTH/ 

1  (FNUZ*VIBPF) 

IF(SR(I) .GT.SMAX)  SMAX=SR(I) 

ER ( I ) =FREQ 


STRE=SR(I) 

IF ( FREQ . GE . FMIN . AND . FREQ . LE . FMAX )  THEN 
CALL  STORE (NLINE, FREQ, STRE) 

ENDIF 

5000  CONTINUE 

DO  6000  I  =  JMAX, JMIN+1, -1 
SREL  =  SP(I)/SMAX 
WRITE(6, 1099)  I , EP (I) , SP ( I) , SREL 
6000  CONTINUE 

DO  7000  I=JMIN, JMAX-1 
SREL=SR ( I ) /SMAX 

WRITE(6, 1109)  I,ER(I) ,SR(I) , SREL 
1109  FORMAT ( '  R  ' , 3X, 14 , 5X, F10 . 4 , 5X, 1PE9 . 3 , 5X, F10 . 6) 
7000  CONTINUE 

WRITE (6, 1119)  TOTINT 
1119  FORMAT ('  TOTINT  =  ' ,F10.6) 

CALL  SORT (NLINE) 

CALL  BINSORT (NLINE) 

CALL  SHAPE (NLINE) 

STO 


Table  l  -  Statistical  Weight  Factors  t'or  Isotopes  of  Water 


Isotope 

Value  of  |r| 

Symmetry 

Relative  Weight 

H,0 

odd 

antisymmetric 

3 

even 

symmetric 

1 

D,0 

odd 

antisymmetric 

3 

even 

symmetric 

6 

HDO 


odd  or  even 


no  identical  nuclei 


Tabie  2  -  Honl-London  Factors  for  CO 
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Table  3  -  Constants  for  H,0  Calculations 


Rotational  Constants  (in  cm"1) 


V1V2V3 

000[9] 

100[10] 

0 1 0[9] 

00 1  [ 1 0] 

020[  10] 

A 

27.88067 

27. 122 17 

31.1 283s 

26.64805 

35.5S6- 

B 

14.521 680 

14.3047, 

14.687569 

14.43130, 

14.8415 

C 

9.277459 

9. 10457 

9.12913, 

9.13816, 

S.9'44s 

AjXlO3 

1.24894 

1.2330 

1.39537 

1.30549 

1.5804 

Ajkx103 

-5.765s 

-5.3874 

-7.605j 

-5.656  j 

- 1 0.48s 

Akx102 

3.25 19g 

3.023o 

5-755,3* 

2.8584 

1 0.99 1 Q 

5jXl04 

5.083a 

4.9987 

5.7870  . 

5.38 1  - 

6. '5, 

5kx103 

1.3007 

1.240s 

3.766, 

1 .326  x 

8-'li 

Vibrational  Constants[8] 

v' 

V* 

i/jlcm'1) 

S®\102O(cm  molecule  1 

Band  T\pe 

100 

000 

3657.053 

48.62 

B 

010 

000 

1594.778 

1038.0 

B 

001 

000 

3755.930 

692.5 

4 

020 

000 

3151.630 

7.537 

B 

Table  4  -  Constants  for  D,0  Calculation 


Rotational  Constants 


V1V2V3 

000[  1 1  ] 

1 00[  1  1] 

0 1 0[ 1 2) 

00 1[  I  1] 

AIcnT1) 

15.3846 

15.1286 

16.633880, 

14.7916 

Blcm'1) 

7.2716 

7.1696 

7.338823! 

7.2266 

Ctcm'1) 

4.8458 

4.7698 

4.789485g 

4.7908 

Vibrational  Constants 

v' 

v"  v0(cm'l)[13] 

S°\1020( 

cm,  molecule )[  1 4] 

Band  Type 

100 

000  2672.0811 

34.0 

B 

010 

000  1178.374 

564. 

B 

001 

000  2787.7182 

486. 

A 
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Table  5  -  Constants  for  HDO  Calculation 


Rotational  Constants 


V1V2V3 

000[15] 

100[1 6] 

01  Op  7] 

00 1  [  1 3] 

020(16] 

Aicm'1) 

23.413842  . 

23.095854 

25.551 

22.3226 

28.25549g 

BIcm'1) 

9.103323 

8.9252 1  s 

9.238 

9.0850 

9.36779- 

Clem'1) 

6.406295 

6.304243o 

6.335 

6.3293 

6.24007-g 

Vibrational  Constants 

v' 

v” 

t'0(cm‘1)[I5] 

SJj’x  1 02O(cm/ molecule ) 

Band  Type 

100 

000 

2723.680 

215.  [18] 

A 

100 

000 

2723.680 

1.73  [18] 

B 

010 

000 

1403.421 

406.2  [19] 

A 

010 

000 

1403.421 

532.8  [19] 

B 

001 

000 

3707.47 

350.3  [20] 

A 

001 

000 

3707.47 

145.3  [20] 

B 

020 

OCo 

2782.012 

24.2  [18] 

A 

020 

000 

2782.012 

11.1  [18] 

B 

*  1,.* 


w-t  tr- '  vri  -  \n-  *  ■  *r.i  yvnn”  yrgrp*  ?XP  n*  n?*n*  t.1"  ivy^r,v»  <  m  *vp  n^w  ^mrrm 


Table  6  -  Constants  for  C02  Calculation 
Rotational  Constants[7] 


V* 

Cv(cm'M 

Blcm'1) 

Dx  107(cm'1) 

H\1013(cm‘ 

00001 

0. 

0.39021817 

1.33204 

0.055 

01101  e 

667.3801 

0.39063825 

1.35133 

01101  f 

667.3801 

0.39125388 

1.35900 

000 11 

2349.1433 

0.38714069 

1.32873 

0.077 

01111  e 

3004.0122 

0.38759172 

1 .34546 

01111  f 

3004.0122 

0.38818943 

1.35522 

10012 

3612.8417 

0.38750237 

1.57314 

2.024 

10011 

3714.7828 

0.38706251 

1.14177 

1.989 

♦For  band  notation,  seep], 

pp.  23-24 

Vibrational  Constants[7] 

u0  (cm-1) 

v' 

v" 

S^.x  1 02O(cm/  molecule) 

Band  Type 

667.380 

01 101 

00001 

82580. 

nu  -  ££ 

2336.632 

01111 

01101 

73700. 

ng  -  nu 

2349.143 

00011 

00001 

955900. 

v*+  0  + 

-u  -  -g 

3612.842 

10012 

0000 1 

10400. 

r>+  v* 

-u  -g 

3714.783 

10011 

00001 

15800. 

-U  -g 
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Table  7  -  Constants  t'or  HC£  Calculations 


Rotational  Constants[24] 


Blcm"1) 

Diem'1) 

H3SC£(v=0) 

10.4400 

0.000530 

H3SC£(v=1) 

10.1381 

0.000526 

D3SC£(  v=0) 

5.392149 

0.0001397 

D35C£(v=1) 

5.278858 

0.0001410 

Vibrational  Constants[24] 


296  K 


Simulation  of  hot  COs  emission  spectrum  at  I730K.,  cold  CO,  absorption  spectrum  at 
296K.,  and  hot  COj  emission  viewed  through  absorbing  cold  CO.. 


